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A Hybrid Numerical Method for 
Turbulent Mixing Layers 


Abstract 


by 


NICHOLAS JACOB GE0RGIAD1S 


A hybrid method has been developed for simulations of compressible turbulent 
mixing layers. Such mixing layers dominate the flows in exhaust systems of modern 
day aircraft and also those of hypersonic vehicles currently under development . The 
method developed here is intended for configurations in which a dominant 
structural feature provides an unsteady mechanism to drive the turbulent 
development in the mixing layer. 

The hybrid method uses a Reynolds-averaged Navier-Stokes (RANS) procedure to 
calculate wall bounded regions entering a mixing section, and a Large Eddy 
Simulation (LES) procedure to calculate the mixing dominated regions. A 
numerical technique was developed to enable the use of the hybrid RANS-LES 
method on stretched, non-Cartesian grids. Closure for the RANS equations was 
obtained using the Cebeci-Smith algebraic turbulence model in conjunction with the 
wall-function approach of Ota and Goldberg. The wall-function approach enabled a 


NASA/TM — 2001-2108 1 1 


XVII 


continuous computational grid from the HANS regions to the LKS region. The LES 
equations were closed using the Smagorinskv suhgrid scale model. 

The h\ brid R A NS- LES method is applied t o a benchmark compressible mixing layer 
experiment. Preliminary two dimensional calculations are used to investigate the 
effects of axial grid density and boundary conditions. Vortex shedding from the base 
region of a splitter plate separating the upstream flows was observed to eventually 
transition to turbulence. The location of the transition, however, was much furt her 
downstream than indicated by experiments. 

Actual LES calculations, performed in three spatial directions, also indicated vortex 
shedding, but the transition to turbulence was found to occur much closer to the 
beginning of the mixing section, which is in agreement with experimental 
observations. These calculations demonstrated that LES simulations must be 
performed in three dimensions. Comparisons of time-averaged axial velocities and 
turbulence intensities indicated reasonable agreement with experimental data. 
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CHAPTER 1 


INTRODUCTION 


The use of computational fluid dynamics (CFD) to assist in the analysis and design 
of aerospace vehicles and their components has substantially increased in recent 
years. For analyzing one particular class of flows, that of aircraft engine exhaust 
nozzles. Reynolds-averaged Navier-Stokes (RANS) codes have been used extensively 
by government organizations (i.e. NASA) and aerospace companies. Exhaust 
nozzles being developed for modern day subsonic commercial aircraft typically have 
multiple streams with a core flow and one or more bypass streams which mix with 
the high energv core flow before exiting the nozzle to lower jet noise while 
maintaining high thrust levels. Similarly in NASA’s High-Speed Research program, 
the engine exhaust systems for the proposed supersonic transport were designed to 
be mixer-ejector nozzles, which entrain secondary air into the exhaust nozzle to mix 
with the core engine stream, again with the goal of simultaneously lowering jet noise 
and maintaining sufficient thrust [91]. Propulsion systems currently under 
development for use on hypersonic and reusable space launch vehicles, such as the 
Turbine-Based Combined-Cycle (TBCC) and Rocket-Based Combined-Cycle 
(RBCC) concepts also employ mixer-ejector ducts. The TBCC concept [30. 104] 
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us< s s a mixer-ejector to integrate a turbine engine with a ramjet into a propulsion 
system with a common nozzle, while the RBUU concept [28. 101. 1 0-'i] integrates a 
rocket with a ramjet, again with a common nozzle. 

The flows in these nozzle systems all have compressible turbulent mixing as the 
dominant flow characteristic. RANS codes used by research and development 
engineers to analyze these nozzles have employed turbulence models to replace the 
unsteady turbulent motion with an effective eddy viscosity. Unfortunately, no 
turbulence model lias been developed to date which is able to accurately represent 
the turbulent motion for such nozzle flows. References [5] and [32] show that the 
“state of the art" turbulence models available in production-use HANS codes have 
major deficiencies in predicting turbulent mixing in nozzle and jet flows involving 
compressibility, high temperatures, and three-dimensionality. 

The known limitations of RANS techniques to calculate complex turbulent flows, 
coupled with continually increasing computing power, have led to interest in more 
sophisticated calculation techniques such as direct numerical simulation (DNS) and 
large eddy simulation (LES). DNS is currently limited by computer hardware to 
very simple flows at low Reynolds numbers, and LES, which directly solves for the 
large turbulent scales and limits empirical modeling to the smallest scales, is 
becoming practical for more complex flows at higher Reynolds numbers. Birch [10] 
and Bradshaw [14] suggest that LES techniques offer the best prospects for 
improving the capability to calculate turbulent flows, particularly for flow regions 
not including wall boundary layers. 
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As a result, ati LES-based technique is an attractive option for calculating the 
mixing dominated regions of nozzle flows. However, applying such an LES technique 
simultaneously to the wall bounded regions that, enter the mixing region (which are 
an important part of multi-stream nozzles that should be calculated accurately) will 
not be practical in the near future. This is because computational resources lar 
greater than those available today would be required to capture the wide range of 
turbulent time and length scales that are important in such a problem. These 
turbulent scales range from very' small eddies in the wall boundary' layers to very 
large eddies in the developing mixing layer. 

While RANS-based methods have major deficiencies in predicting compressible 
mixing layers and inherently are not formulated for calculation of unsteady 
turbulent flows, they' have been shown to predict the mean flow behavior of wall 
bounded regions quite well, particularly in the absence of adverse pressure 
gradients. As a result, it would be desirable to combine a R ANS-based technique for 
the wall boundary layers upstream of the mixing region with an LES-based 
technique for the downstream unsteady, turbulent mixing region. The development 
of such a hybrid RANS/LES approach is the subject of this work. 

In the rest of this chapter, a survey' of methods currently used for the computational 
modeling of nozzle flows dominated by compressible turbulent mixing and new 
techniques under development are presented. First, RANS methods are discussed in 
section 1.1 with emphasis placed on turbulence modeling, which is most frequently 
considered the limiting factor of such CFD simulations. Next, a discussion of LES 
techniques, as applied to nozzle and jet flow fields, is presented in section 1.2. 
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Section l.:{ discusses hybrid methods, which have recently been proposed for flows 
in which wall boundary layer regions may be adequately calculated with a HANS 
method, while unsteady regions with large scale recirculations are calculated with an 
EES method. An overview of the hybrid method developed in this work is presented 
in section 1.4. Finally, an outline of this dissertation is provided in section l.b. 

1.1 RANS Methods 

The RANS codes used to simulate the nozzle flows discussed at the beginning of 
this chapter are largely general-purpose flow solvers, capable of handling a variety of 
turbulent flows extending beyond only nozzle and jet problems. State-of-the-art 
turbulence models, such as Reynolds-stress closures, have found their way into codes 
intended for basic research and simple flow problems, but the general-purpose codes 
used for nozzle simulations and other realistic configurations usually employ 
computationally cheaper eddy viscosity models. Eddy viscosity models use the 
Boussinesq approximation to calculate the Reynolds stress as the product of an 
eddy viscosity and the rate-of-strain tensor. These models vary in complexity from 
algebraic (also termed zero-equation) formulations to more complex formulations 
which solve additional transport equations (normally one or two partial differential 
equations). 

Algebraic models typically use a form of PrandtEs mixing length hypothesis to 
calculate the eddy (or turbulent) viscosity, but are usually optimized for a single 
flow and are not accurate for a wide range of flows. One of the first widely used 
algebraic models was developed by Cebeci and Smith [18,19]. Another widely used 
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algebraic model is that due 1o Baldwin and Lomax [1]. Both the ( ebeci-Smilh and 
Baldwin- Lomax models are formulated only for the calculation of wall boundary 
layers and have been used successfully for calculations of flows without adverse 
pressure gradients or separation regions. 

One-equation models offer, to some degree, more generality than algebraic models, 
because they solve for a quantity such as the turbulent kinetic energy or the 
turbulent viscosity. However, they frequently have the limitation of not solving a 
transport equation for the turbulent length scale. Two of the more widely used 
one-equation models are the Baldwin-Barth [3] and Spalart-Allmaras [98] models. 
While one-equation models have shown promise for turbulent wall bounded flows at 
a lower computational cost than algebraic models, they have not been shown to be 
accurate for flows with significant turbulent mixing, which is of piimaiy impoitance 
to the nozzle problems discussed here. 

Two-equation turbulence models usually solve one equation for the turbulent kinetic 
energy' and a second equation from which the turbulent length scale can be 
obtained. The most popular of these are k-e models, that solve one transport 
equation for the turbulent kinetic energy, k, and the second for the rate of turbulent 
kinetic energy dissipation, e. Probably the one k-e turbulence model that has served 
as the basis of most other k-e models (and itself is still in wide use) is that due to 
Jones and Launder [47]. Also referred to as the “standard k-e model, the 
Jones-Launder model has one form that may be directly integrated down to solid 
surfaces including the laminar sublayer, through the use of damping terms, and a 
second form that requires the use of a wall function to bridge the gap from the fully 
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turbulent- region of a boundary layer to t he laminar sublayer very near t he wall. 

1 he former is referred to as the "low Reynolds number** form while the latter is 
referred to as the "high Reynolds number*' or "wall-function** form. Both have the 
same formulation away from walls and in mixing layers. 

Although k-t models are generally more accurate than the simpler algebraic or 
one-equation models, they have several limitations. For wall bounded flows, k-e 
models are particularly deficient in regions of strong adverse pressure gradients and 
in separation zones (Rodi and Scheuerer [84]). For high-speed flows involving 
mixing of multiple streams, such as those in the nozzle and jet flows described 
previously, standard k -c models deviate substantially from experimental data for 
highly compressible shear layers and for round jets. 

Several new k-e: models have been developed beyond the standard Jones- Launder 
model to address pressure gradient and separated boundary layer problems, but 
have not demonstrated wide spread improvement over the Jones- Launder 
formulation. To address the compressibility issue, several modifications have been 
proposed. Most of these modify the equation for dissipation rate. e. so as to account 
for the experimentally observed reduction in turbulent kinetic energy production 
with the relative speed of streams (frequently referred to as the convective Mach 
number) that form compressible shear layers. The mostly widely used are those of 
Sarkar [88,89] and Zeman [ITS]. 

Experimental data have also indicated that round (or axisymmetric) jets mix less 
rapidly than planar jets. Because the empirical coefficients of most k-c models. 
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including Jones- Launder, have been calibrated against incompressible planar shear 
layers, corrections have been developed to improve the capability of k-< models to 
calculate round jets. The most widely used of these is that due to Pope [i9], which 
like the compressibility corrections of Sarkar and Zeman. modifies the dissipation 
rate equation. Essentially, the Pope correct ion increases the dissipation rate of 
turbulent kinetic energy in the presence of vortex stretching, which is characteristic 
of round jets. Finally, one recently developed k-r model, due to Tines and 
Tam [102], substantially modified the empirical coefficients from those of the 
standard k-c model, and included the Sarkar compressibility correction and Pope 
round jet correction intended for turbulent jets with flow conditions similar to those 
found in the exhaust nozzles of high speed aircraft. The Tam-T flies model, with its 
substantially modified coefficients, has only been calibrated for mixing layers using 
calculations beginning downstream of any nozzle surfaces. As a result, some other 
method is required to accurately calculate the wall boundary layer regions upst ream 
of the mixing layer. 

Despite some of the aforementioned limitations of k-c models and the necessity of 
empirically based corrections to handle specific flow- complexities, k-c models are 
still the most frequently used models for calculating nozzle and jet flow' fields. This 
is because more complex models, such as Reynolds stress closures, have not 
demonstrated significant accuracy improvements that would warrant their much 
more computationally expensive use. As a result, several validation studies have 
been conducted to determine the accuracy of k-t models for specific classes of nozzle 
and jet problems. Some recent examples of such validation studies are provided in 
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references [5,0.2*2*31, /5. 106]. Nearly all of the flows investigated have operating 
conditions (velocities and temperatures) similar to those of real configurations, but 
the geometries were simpler to allow comprehensive studies of different models and 
corrections to bo conducted. 

Two other popular two-equation models will bo briefly discussed here. The k-w 
model is another two-equation turbulence model which is similar in form to the k-r 
model. The second transport equation of this model solves for the dissipation rate 
per unit turbulent kinetic energy. The most widely used k-u: in current use is that 
due to Wilcox [113]. The k-u: model has demonstrated improved capability to 
handle adverse pressure gradient boundary layers relative to the k-< model. 
However, it has been shown to be worse than k-( in predicting free shear layer 
mixing, and demonstrates significant sensitivity to freest ream turbulence levels. 

Wi lcox [11-1] extended his original model in a new formulation to address these 
limitations. The last two-equation model discussed here is the Shear Stress 
Transport Model (SST) of Menter [66.67]. Mentor's model employs a k-u: 
formulation in regions near solid boundaries and a k-c formulation away from walls. 
This model has become popular in recent years because it used the k-u; and k-e 
formulations in the regions where each performs the best. 

The most recent adaptation of a Reynolds-averaged Navier-Stokes approach has 
been in the emerging field of computational aeroacoustics to address the problem of 
developing quieter aircraft engine exhaust nozzles. One procedure is to use the 
solution from a turbulent CFD simulation of a jet flow field (using a k-t model) as 
input into an acoustic prediction method. The jet noise field is then determined by 
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integration of the sound propagation equations. References [50. >1] present some 
recent aeroacoustic calculations of convergent-divergent nozzles having flow 
characteristics again representative of realistic high speed nozzles. In computational 
aeroacoustics. which has its own set of limitations related to the acoustic modeling 
methods, the input is also limited by the quality of the aerodynamic flow field. As 
with the nozzle and jet aerodynamic prediction methods, significant work is being 
performed to improve the accuracy of acoustic predict ion techniques [06]. 

1.2 LES Methods 

Although LES techniques have been used in selective applicat ions, such as 
meteorology and atmospheric sciences [S3] for several years, only recently has LES 
been considered as a potential tool for nozzle and jet flow problems. Much of the 
foundation of LES was established in the meteorology field by Smagoi insky [9-1], 
Lilly [58], and Deardorff [25] near the end of the 1960’s. Reference [69] mentions 
that after some interest in LES by the United States engineering community in the 
early 1970's, the 1980’s saw only limited development and application. Recently, 
however, the realization of the shortcomings of RANS, and improvement in 
computer speeds have sparked new interest in LES for engineering applications. 

Several recent U.S. aeronautics research and development programs have sponsored 
research in developing and applying LES based techniques for analyzing jet flows. 
Of particular interest to the these programs is the potential of LES as a 
computational aeroacoustic technique for enabling accurate predictions of jet noise. 
In reference [65], Mankbadi et al. discusses that large scale turbulent stiuctuies are 
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believed to be the dominant noise producing mechanisms in free jet flows. Choi et 
al. [20] have use an LES technique to calculate the aerodynamic field of a Mach 1.1 
axisymmetric jet and then employed a Kirchoff method, detailed in references [00] 
and [02], to calculate the noise field. Since LES reserves modeling for the smallest 
turbulent length and time scales while directly solving for the large scales, it is 
thought to be appropriate for applicat ion to certain classes of jet noise problems. In 
particular, success is being realized for jet problems in which the flow regions of 
interest are far removed from solid boundaries, such as relatively simple single flow 
nozzles exiting into ambient air. LES methods are not vet prepared to handle more 
complex configurations, like those of the multiple stream nozzles discussed 
previously. Such configurations, if they contain turbulent boundary layers along 
solid surfaces, have a large range of turbulent scales that are important to the 
overall flow problem, and at this point , t he cut-off scale at which direct solution 
ends and modeling begins is one of the major unresolved issues of LES. 

A brief review of the fundamental steps of LES is as follows, more details of the 
process will be given in chapter 2. A filtering process removes the smallest spatial 
scales of the Navier-Stokes equations. A set of equations is produced that represents 
the spatial and time evolution of the larger turbulent scales, but contains a subgrid 
scale tensor to account for the unresolved smallest scales that were removed by the 
filtering process. Next, the subgrid scale stress is replaced by a model 
(appropriately termed a subgrid scale model), which is analogous to the eddy 
viscosity models described in the last section for Reynolds-averaged Navier-Stokes 
codes. I sing the subgrid scale model to handle the small scale structures ranging 
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from the Kolmogorov scale up to the cut-off scale, the large scale equations are then 
solved using a time-accurate algorithm. 

The features of the subgrid scale model has been one of the most heavily researched 
aspects of LES. The most commonly used subgrid scale model is the Smagorinsky 
eddy viscosity model [91]. One of the primary reasons for t he popularity ol the 
Smagornisky model is that it provides sufficient diffusion and dissipation to keep an 
LES computation stable. Because of the importance of the subgrid scale model in 
allowing for the correct transfer of turbulent energy between the large scales (which 
are directly solved) and the modeled subgrid scales, several modifications to the 
Smagorinsky model have been proposed. The most popular of these are known as 
dynamic eddy-viscosity models [26.33,77] which replace the model constant of the 
Smagorinsky model with a coefficient that is allowed to vary both spatially and 
temporally as a calculation progresses. The major advantage of the dynamic 
eddy-viscosity models over the Smagorinsky model is the improved capability of the 
dynamic eddy-viscositv models to provide the correct turbulent kinetic energy 
dissipation. As with the turbulence models for Reynolds-averaged Navier-Stokes 
codes, however, no single subgrid-scale model has been found to work well for all 
flows. Nelson [70] and Vreman [107,108] have investigated the effects of subgrid-scale 
model selection on LES calculations of compressible planar shear layers. 

There are other issues related to LES that will require significant work to make LES 
an appropriate tool for engineering applications, and in particular foi the nozzle and 
jet flows that are the subject of this report. Since LES directly solves for the largest 
scales of turbulent motion, a highly accurate solver is necessary' for the 
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time-marching procedure. References [39.55.68] discuss higher order schemes as 
applied to LES. Compact finite difference schemes, introduced by Lele in 
reference [55]. are particularly promising for LES calculations because they are 
available' to accurately resolve a greater range of scales t han other finite difference 
schemes. Several LES algorithms are using (Mae( ormack-type) explicit schemes, 
which lend themselves well to calculations distributed over parallel processors (see 
references [10] and [93]). Boundary conditions provide some of the most challenging 
modeling uncertainties for LES. Many LES simulations of simple geometry' 
benchmark cases, such as fully-developed channel flow, have used periodic boundary 
conditions for the computational inflow and exit. For jet flows, however, such a 
boundary condition cannot be used, because of the fundamentally different nature 
of the inflow and outflow stations. 

Nearly all LES simulations of jet and mixing layer flows performed to date have 
placed the inflow of the computational domain downstream of any wall bounded 
regions and have either ignored the upstream boundary layer effects or used some 
approximation to initialize the turbulent mixing layer. Several authors, such as 
Ragab [81,82] and Hedges [38] have imposed hyperbolic tangent mean velocity 
profiles at the plane which represents the end of the wall boundary layer regions and 
the beginning of the mixing region. Ragab then used the results of a linear stability 
analysis to generate perturbations about the mean velocity profile located at t lie 
mixing plane. Hedges added small amplitude perturbations to the vertical velocity 
component in simulations of heated jet flows. 


NASA/TM — 2001-210811 


12 



The difficulty in using an artificially generated inflow, such as that assuming a 
hyperbolic tangent velocity profile, is that the characteristics of the upstieam 
turbulent boundary layers, including velocity, temperature, and turbulence profiles, 
are not accurately represented. This is a. significant deficiency since the state of the 
incoming boundary layers have been shown to have significant effects on the 
development of turbulent mixing layers in the experiments conducted by 
Bradshaw [13], Browand and Latigo [15], and Hussain arid Clark [44]. 

Recently, Li et al. [57] proposed a method in which parallel calculations of the 
upstream boundary layers are used to generate time-varying inflow conditions foi a 
spatially developing mixing layer. This method offers a. promising technique for 
reducing the computational cost relative to performing a complete LES calculation 
containing both the wall-bounded and mixing regions, but itself is probably too 
expensive to use in the near future for high speed, high Reynolds number cases. The 
hybrid RANS-LES method developed in this work is proposed as an alternative 
computational technique to performing LES calculations everywhere in the 
computational domain, that includes the mean flow characteristics of the incoming 
boundarv layers and is also feasible when considering foreseeable computational 
resources. 


1.3 Hybrid RANS-LES Methods 

The realization that LES calculations of flows in aerospace and industrial 
applications at realistic Reynolds numbers will not be possible for some time has led 
to interest into the development of hybrid techniques. The objective of a hybrid 
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me) hod is to retain the essential feat ures of the LES method, but to employ a 
computationally cheaper RANS method in regions where it is appropriate. As a 
result, nearly all hybrid methods proposed to date apply a RANS approach to 
attached wall boundary layer regions and an LES approach to regions of large scale 
separation. The work detailed in this dissertation represents the first hybrid method 
development for applicat ion to compressible mixing layers. 

The most widely publicized hybrid method to date is the Detached Eddy Simulation 
(DES) method of Spalart [96,97,99]. In the DES method, the wall bounded regions 
are calculated using RANS with the Spalart- A llmaras [98] one equation turbulence 
model. Constant inescu and Squires [23] have applied Spalart \s DES method to 
turbulent flow over a sphere, which is an appropriate geometry for the method due 
to the large scale separation in the wake of the sphere. 

Speziale [100] suggested an approach that allows for computations varying from 
RANS in the coarse grid limit, through LES, and finally to DNS in the very fine grid 
limit. A Reynolds Stress model is used to close the turbulent stresses in the RANS 
limit, and provides the basis for a subgrid model necessary in LES simulations. 
Batten et. al. [7] also propose a hybrid model that employs a Reynolds-Stress model 
to close the RANS and LES equations. Lastly. Arunajatesan et al. [1] have applied 
a hybrid RANS-LES method to cavity flovvfields. Their approach employed a 
two-equation k-kl turbulence mode] to close the RANS equations and a one-equation 
model solving for the filtered subgrid kinetic energy to close the LES equations. 
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1.4 Overview and Objectives of the Current Work 


The hybrid RANS-LES method presented in t his dissertation is developed for 
application to configurations such as the mixing layer shown in figure 1.1. This 
relatively simple configuration is representative of the more complex nozzle 
geometries discussed at the beginning of this chapter, in that two wall bounded 
regions provide isolated flows to a single region where compressible mixing is the 
primary flow characteristic. Development of the hybrid method and assessment of 
the method for a benchmark compressible mixing layer configuration are the focus 
of the dissertation. 

The hybrid method employs a BANS approach to provide the mean flow 
characteristics of the wall boundary layers entering the mixing region. The 
downstream mixing laver is then calculated using LES. The method developed heie 
is intended for those nozzle and mixing layer problems in which a dominant 
geometric feature, such as the base region of a nozzle or splitter plate separating the 
upstream flows, will provide an unsteady mechanism to drive the tuibulent 
development in the mixing layer which will dominate unsteady effects from the 
incoming turbulent wall boundary layers. Although the upstream RANS approach 
does not provide any unsteady turbulent information to the mixing layer, the mean 
flow momentum and thermal boundary layer effects can be provided. 

1.5 Outline of the Dissertation 

The equation sets which must be solved for the RANS and LES regions are derived 
from the general form of the compressible Navier-Stokes equations in chapter 2. The 
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Wall Bounded Regions - RANS 


Free Shear Layer Region - LES 


Figure 1.1: Schematic of mixing layer demonstrating the hybrid RANS/ LES ap- 
proach. 


turbulence model used to close the RANS equations and the subgrid scale model 
needed to close the LES equations are also presented. Chapter 3 provides details of 
the numerical procedure used to perform the hybrid RANS-LES computations. 
Validation of the RANS method for a series of flat plate boundary layer cases is 
presented in chapter 4. These cases include an incompressible laminar boundary 
layer, an incompressible turbulent boundary layer, and finally two supersonic 
boundary layers which have the same flow conditions as the two streams entering 
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the compressible mixing layer that is the focus of this work. 1 he lnbtid HANS EES 
method is investigated for the compressible mixing layer in chapters •) and (>. 
Although true LES simulations require calculations in three spatial directions, 
two-dimensional calculations are examined first in chapter 5 to investigate effects of 
grid resolution and boundary conditions. Three dimensional calculations are 
examined next in chapter (> with emphasis placed on comparison to experimental 
data, and to parametric st udies of subgrid scale model and grid resolut ion effects. 
The importance of using three dimensional calculations to capture the initial 
development of the turbulent mixing layer is also investigated. C ontlusions of the 
hybrid RANS-LES method development and benchmark computations as well as 
recommendations for further research are presented in chapter 7. 
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CHAPTER 2 


FORMULATION OF THE RANS AND LES EQUATIONS 

In this chapter, the equations that are used to solve the flow in the HANS and LES 
regions are derived. The general form of the Navier-Stokes equations, written in 
tensor notation, is presented first. Next, the RANS equations are derived using the 
mass-weighted form of the Reynolds-averaging process and the LES equations are 
derived using a mass-weighted spatial filtering process. Turbulence modeling foi the 
RANS and LES equations is presented in the last section of the chapter. 

2.1 Navier-Stokes Equations of Motion 

The Navier-Stokes equations represent the time-dependent, three-dimensional 
motion of a fluid. They consist of expressions for the conservation of mass, 
momentum, and energy. 

The expression for conservation of mass, or the continuity equation is written in 
tensor form as: 


% + (P u *‘) = 0 

dt Ox, 


( 2 . 1 ) 
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Conservat ion of momentum is written: 


0 () OP Or, 

(/-»«/) + TJJ- ~ + 


ot 


Ox, Ox, 


{■):, 


:l) 


Conservation of energy is expressed as follows: 


[Uj ( E ' + I ’ )) -J7 J ( “' T *' ) - wp 

Hen', the variable E, represents the total energy (internal energy plus kinetic 
energy) per unit volume: 


(2.4) 


E, — pc + - pu,u , (2.4) 

The equation of state for an ideal gas is used to relate the pressure, temperature, 
and density through: 


P = pHT 


(2.5) 


For the viscous stresses r,j. it is assumed that the fluid is a Newtonian fluid, and as 


a result , the viscous stress is proportional to the rate of strain. This is written: 


Tjj — 2pSij + A— 


du j 
'Ox/ 


where the rate of strain tensor S,j is: 


r. _ I f9ui , e)u j 
" ‘ J 2 \ Ox, Ox, 


( 2 . 6 ) 


(2.7) 
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Using Stokes's assumption that the thermodynamic and mechanical piessuies aie 
the same for a fluid undergoing and expansion or compression: 



( 2 . 8 ) 


Equation (2.6) can then he rewritten as: 


Tij = 


2 Ch, J r 

d.Tj 


3 


(2.9) 


To calculate the viscosity, the Sutherland model is used which assumes that the 
viscosity is only a function of temperature for a gas: 


ft = 


U.7- 

C 2 + 7 


( 2 . 10 ) 


For air. the values of the constants (\ and C 2 are (in SI units): 


C’i = 1.458 x 10 
C 2 = 110.4 A 


-6 


kg 


rn ■ $ ■ I\ 2 


( 2 . 11 ) 


The heat flux g, is obtained from Fourier's law: 


where k is the thermal conductivity. It is assumed that the fluid is thermally 
perfect, such that the internal energy and enthalpy are only functions of the 
temperature, and it is also assumed that the fluid is caloric-ally perfect such that the 
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specific heats ( V and C P are constants. As a result, the interna] energy and 
enthalpy ('an he written as: 


< =('\T 

h = C P T 


( 2 . 13 ) 


Assuming tfiat the air is of constant composition and does not undergo any 
chemical reaction, the thermal conductivity is only a function of temperature. Using 
the specification of constant specific heats, the following expression is obtained for 
the t hermal conductivity as a function of the constant pressure specific heat, 

Prandtl number, and the viscosity defined in equation (2.10): 


2.2 Mass- Weighted RANS Equations 


In the classical form of Reynolds averaging, the time dependent form of the 
Navier-Stokes equations given by (2.1) through (2.3) are averaged over a period of 
time that is much larger than the period of turbulent fluctuations. Each of the 
dependent variables appearing in these equations is replaced bv the sum of mean 
and fluctuating components. As an example, the velocity would be given by: 


Ui = U{ + u- 


(2.15) 


where the time averaged velocity tq is given by: 


1 /H-r 
Ui = - / u 
T Jt 


(It 


(2.16) 
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For the current work, where fluctuations in density are important, a mass (or 


density) weighting is employed in the averaging process, which will make the final 
form of the RAXS equations much more convenient to work with. The dependent 
variables are again broken into mean and fluctuating components: 


where the time averaged (using mass weighting) velocity u, is given bv: 

1 f' +T 

iii — — pu,dt (-• 

pr Jt 

This mass-weighted Reynolds averaging process is frequently referred to as Favre 
averaging, and in general, the Favre average of any variable f is defined b\ . 



(2.19) 


Before averaging the Navier-Stokes equations (2.1 - 2.3) the flow variables are 
separated into mean and fluctuating components: 
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p = p + p 
p = T + p' 

«i = »/ + (/•' 

<- = <+<" (2.20) 

h = /> + h" 

T=f+ T" 

<i> = <ii + <i, 

Note that the density, pressure, and heat flux are not decomposed using mass 
weighted variables. Starting with the continuity expression (2.1). a time averaging is 
performed to obtain: 


dp 

dt 


+ (P u i) — 0 
o-r, 


1 1 sing the definit ion of mass weighted variables, this can be rewritten as: 


( 2 . 21 ) 


dp d 

- + s -( ?a ,)=° 

Working next with the momentum equation (2.2). a time average of the entire 
equation is performed, resulting in: 


In equation ( 
variables as: 


0 () dP Ot j j 

t - {pun + — (pupij ) = - — + 


dt 


c)x , 


Ox, OXj 


(2.23) 


23). the time dependent term is rewritten in terms of mass averaged 
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J, (n«i i = A (?«.-) 


The convective term is expanded as: 


77— (/>"/",) = -7— (p(«» + "D ("./ + "” )) 

O.V j ' OXj v 

= Jy (w) + ^ + d7 t 

* ,/ 7 J 

The first term in the expansion of (2.25) is rewritten as: 


77 - {(riiiUj) = 77 - (/atift/j) (2.26) 

d.Vj v rAr; 

Over the period of the time averaging, the the mass weighted variables denoted with 
a hat are constant. As a result: 


« (^0 


(2.27) 


The last term in equation (2.25) is the turbulent or Reynolds stress, and is the term 
from the momentum equation that is replaced with a turbulence model. 


0 / , dr? 

{pv ' u ' ;) = “a 7 


(2.28) 


Replacing the terms in equation (2.23) with those in equations (2.24) through 
(2.28), the resulting RANS momentum equation written in terms of mass averaged 


variables is: 
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To obtain the time averaged energy equation written in mass averaged variables, the 
total energy terms in equation (2.3) are replaced with the expanded energy 
expression shown in (2.4) and then the entire equation is time averaged: 


(pc + tpu.u,) + ~ (f Ij {pt + \pupp + P)) = 

Working first with the time dependent term, this term i 
plus fluctuating velocities as: 


8 i — i i)q j 


dx : , 


dx 


(2.30) 


s expanded in terms of mean 


£ (p< + ^ (pc) + £ (\p(ui + «") (u, + <)) 


(2.31 


The kinetic energy term is further expanded as: 


J t (Jp (“>• + «?) («i + «n) = (hpUiU.) + ~ (\pipu") + ^ (kp«) (2.32) 


The the first term on the right side of equation (2.32) is rewritten as: 


d_ 

dt 



d_ 

dt 


{\pupp) 


(2.33) 


Over the period of the time averaging, the mass weighted variables are treated as 
constants and as a result, the second term on the right side of equation (2.32) is 
zero: 
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I (!/>*<■'?) = 0 


(2.34) 


The last term is the turbulent kinetic energy and is written as: 



The resulting time dependent term is: 


d_ 

m 


(pc + \plliUi) 



Working next with the convective term, h = c + P/p is used to obtain: 


77 [uj (pc -f \pu,u t + P)) = — (uj {pit + ipu,;/,)) 

OXj v 

Expanding this expression in terms of mean plus fluctuating velocities as: 

(«; (ph + \ pUiUi )) = (p («j + «") (*> + / '") ) 

+ -p- (|p («j + ) («i + "") ("< + U D) 

Working first with the enthalpy expression: 

A (p (% + u") (A + /*") ) = ~ (pujh) + ^ (pujk") 

UlIj J J 


(2.35) 


(2.36) 


(2.37) 


(2.38) 


(2.39) 
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The second and third terms on the right side of equation (2.39) are zero by the 
same argument as that for (2.31). I lie last term contains the turbulent heat flux 
and is rewritten as: 


<) 




()x , 


(2.40) 


I his turbulent heat flux must be modeled, and is done so in a manner similar to 
that for the Reynolds stress appearing in the momentum equation (2.28). 
Working next with t he second term on the right side of equation (2.38): 


() 


("./ + u j) i u > + u '') {<<' + "")) = — 


()Xj \ f ■ •" ' ■ ■ • ■ "/ ()x 


+ 4 ~ (K^f) + (>X“f) 


r) 


dxj v ‘ w ■ Ox, 


Fhe second and fourth terms on the right side of equation (2.41) are zero over the 
period of the time average. Fhe third term includes the turbulent kinetic energy as 
defined in equation (2.35) and is rewritten: 



(2.42) 


The fifth term includes the turbulent Reynolds stress as defined in equation (2.28) 
and is rewritten: 


_d_ 

dxj 


(pUjitjuf) 



(2.43) 
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Tlie resulting convective term is: 


c) 


2.44) 


("> (l ’ h + W*4) = (?»>) + 

+£- (^">'^1 + ^ <*"«) + £-. (w«) 

The viscous dissipation term on the right side of equation (2.30) is examined next 
After expansion, this term becomes: 


d 

d- r j 


( UiTij ) = —— {{ill + ■«?) ( T P + '/,)) 


fix 

d 


c) 


(2.45) 


-!!;«**> +s; w * ) 


The heat flux term is left as it appears in equation (2.30). Rewriting the energy 
equation (2.30) with the terms we obtained in equations (2.36). (2.44). and (2.45) 
results in: 


(pf + \piiiiii + pfy + {piijh + i pujUiUj + piijk) 

= TT- {ilJij + il^fj + UiTij) - -7T- (qj + qj + \p u " u " u ") 

ox j oxj 

In order to simplify this expression to the form that will be used in the 
computations, the terms involving the turbulent kinetic energy and those with the 
fluctuating velocity u" are assumed to be small compared to the other terms in the 
energy equation. The enthalpy in the convective term is first expanded using 
ph = pe + 7. Equation (2.4) is then used to rewrite the sum of the internal and 
mean flow kinetic energy, appearing in both the time-dependent and convective 
terms, as the total energy: 


2.46) 
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E, - />< + 


( 2 , 1 7 ) 


The resulting final form of the RANS energy equation is: 


d_ 

Wt 


(A) + 


— (ujE t + Uj P) 

(JJ j 


0 

dir, 


{iljTij + f/,T,') 


d_ 

c).r, 


fe + 9j) 


2.3 Spatially Filtered LES Equations 


(2.48) 


To derive the LES equations used in this work, the time dependent form of the 
Navier-Stokes equations given in equations (2.1) through (2.3) is used as the 
starting point. Instead of time averaging these equations, however, an approach 
similar to the work of Ragab and Sheen [81.82] and Erlebacher et al [27] is used 
that will filter out small scale fluctuations, and only retain scales that are large 
enough to be resolved by a particular computational scheme and the computational 
mesh. The filtering operation is defined on any variable / bv the expression: 


.f(xJ)= f D G(x-Z.&)f(Z.1)d 3 Z (2.49) 

In equation (2.49), G is the filter function, D is the flow domain, and A is the filter 
width. The filter width A is usually taken to be the grid spacing, and is the 
approach taken here. Note that the overbar used in equation (2.49) indicates a 
filtered variable. This is in contrast to the previous use of an overbar to indicate a 
time averaged quantity in the previous section. 
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As discussed in [71], the exact form of the filter function is not typically known. 
However, the filter function must satisfy: 


/ G {x - £. A) = 1 (2.-)0) 

Jd 

In addition, the form of the function G allows the operations of filtering and 
differentiation to commute such that: 


dt ()t 


&L = K (2.52 

Ox j dx 3 

In large eddy' simulations of compressible flows, it is common to use Favre-filtei mg 
which is defined as: 


1-4 

P 


(2.53) 


w 


.-here a quantity / is decomposed into resolved and unresolved (also referred to as 
sub-grid scale) components as: 


/ = / + /' ( 2 - 54 ) 

Equations (2.53) and (2.19) are very similar in appearance, but they refer to very 
different operations. Both operations employ mass (density) weighting, but equation 
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(2. 19) defines a time averaged quantity, while equation (2.53) defines a spatially 
filtered quantity. As was the case for the RANS equations, the density, pressure, 
and heat flux terms are not decomposed using mass-weighting. Again note that the 
overline represented a time averaging process in the previous section, but it will 
refer to a spatial filtering operation in the current sect ion. In addition, Favre 
filtering differs from Favre time averaging in that: 

/ ± f (2.55) 

and 


r / o 


(2.56) 


Applying the filtering operation to the continuity expression (2.1) results in: 


fjp 

c)t 


+ (pUi) — 0 

Or, 


[Z.oi 


This is rewritten in terms of Favre-filtered variables as: 



3 



Filtering the momentum equation yields: 


(2.58) 




i 2L + f 2hi 

dxi dxj 


(2.59) 
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Rewriting this in terms of Favre-filtered variables results in: 


d () dP (h,, 


The convective term is further expanded: 


2.60 ) 


~r~~~ IpUiUj) = ~~~~ ( p(iii + u’i) {uj + 
aXj ' ox j 

() 


(2.61 


= 7j— [p (ti.-iij + iipij + n'p'ij + »'"')] 

Ol j 

This is reorganized to separate the resolved convective term from the unresolved 
(subgrid-scale) terms as follows: 


{puptj) =~ pUh'ij) + - . [p («;«/ - + 


^-[p («,»' - »/<«,)] I (pu 1 ^) 


(2.62) 


The right side of equation (2.62) contains the resolved convective term. Leonard 
stress, cross stress, and subgrid scale Reynolds stress, respectively. This is often 
written in more compact form as: 


(pUiUj) = TT P + -7T- \p( u i u J - )] (2.63) 

OXj (JXj • OX J 

The cross stress and subgrid stress are frequently modeled together, while the 
magnitude of the Leonard stress is on the order of the truncation error when 
employing most finite difference schemes, and as a result, is implicitly 
represented [114]. Other authors such as Choi et. al [20] neglect the cross stress and 
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Leonard stress contributions entirely. For the method presented in this work, it is 
assumed that the Leonard, cross stress, and subgrid-scale Reynolds stress can be 
modeled together as the subgrid-stress term, as is doin' in reference [17] or: 


' u 


P ( l V«j - I'liUj) + p (««■«'• - ll'p'lj) + pu'pi'j 


(2.64: 


= plitjUj - tijtij) 

1 lie subgrid scale Reynolds stress is replaced with a model. In contrast to the 
Reynolds stress resulting from the time averaging process appearing in equation 
(2.28), in which all turbulent motion is replaced bv a model, the subgrid model 
employed here only represents t he net effects of turbulence motion t hat is not 
resolved by the computational method. More details of both the RANS and LES 
turbulence modeling will be provided in a later section. 


The final form of the filtered momentum expression that will be used for the LES 
calculations is: 




OP dfi, QtI? s 

— j- — — + , H 

<d.r, dxj dxj 


(2.65) 


Many different forms of the energy equation are used by researchers investigating 
compressible flows with LES. Piomelli provides a comprehensive summary of several 
forms ol the energy equation in reference [78]. The form used in this work solves for 
the total energy, as defined in equation (2.3). 


The Favre- filtered energy expression has the form: 


NASA/TM— 2001-210811 


34 



(pt- 4 -Jni,u t ) 4 (ptijh + ^ m j u ' lli ) 


_0_ 

dx, 






( 2 . 66 ) 


The kinetic energy part of the time dependent term is rewritten analogously to the 
convective term of the momentum equation: 


\pu,Ui = \piiiUi 4 \p{uiUi - iiiii-i ) 

(2.67) 

— ^pttjU, 4 k 

where k is the unresolved kinetic energy term. 

The expression involving enthalpy in equation (2.66) can also he lewiitten as the 
sum of resolved and unresolved components: 


pUjh = piijh + p (iijh — Ujh'j (2.6b) 

The unresolved term on the right side of equal ion (2.68) is the subgrid scale heat 
flux, so that equation (2.68) becomes: 


piijh = piijh 4 q'‘ gs 


(2.69) 


The kinetic energy part of the convective term in equation (2.66) is expanded as 


follows: 


kpUjltjUi =±p{UjUiUi 4 2 UjU t u' t + UjU'pi\ + u'jUjUi 
+ 2 u'jUiU'i + ) 


(2.70) 
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Rewriting this in terms of resolved and unresolved components results in the 
following expression: 

^pUjUiUj = 3 p{ujUiiii)+^p(UjUjii , - iijUjit; + : 2t'/, ('/,(/' + u'jUiii, + 

_____ _ ' ( 2 . 
+11.1*1',"', + 

Of unresolved terms in equation (2.71). it is common practice to ignore all terms 
except for the last two term, as demonstrated in reference [82]. Following Hagai), 
the next-to-last term is rewritten as: 


piiju'u • = iijk 

and the last term is rewritten as: 


pUjUiUj 


U;T, 


(2.73) 


The viscous dissipation term on the right side of equation (2.66) is rewritten as the 
sum of resolved and unresolved components as follows: 


d 


rr—( UjTij) = — 
dxj dxj 


d d _ 

( n,T, j ) T — — ( UjTij ** i 1 ij ) 


(2.74) 


It is common practice to assume that only the resolved dissipation term on the right 
side of equation (2.74) is important, although Erlebacher [27] and Piomelli [78] 
indicate that the unresolved contributions should not be ignored for highly 
compressible flows. 
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The heat, flux term Is left unchanged, as was done for the HANS energy equation in 
the preceding section. Rewriting the filtered energy equation (2.66) with the terms 
derived in equations (2.67. 2.69. 2.72. 2.73. and 2.74), the following form is obtained: 


c)l 


{()< + + k) + — (pujh T -T iij k ) 

O.l j 

= {“iTij + -UiT-f' ) - — (<ij + qf ) 


2. / •) I 


The final form of the LES energy equation is obtained using a procedure analogous 


to that used to obtain the final HANS energy equation (2.48). First, the terms 
involving the unresolved turbulent kinetic energy are assumed to be small compared 
to the other terms in the filtered energy equation, and the enthalpy term is recast as 
j)\, — p<' T. Then, the resolved total energy is defined to be the sum of the 
resolved internal energy and the resolved kinetic energy: 


Et — p 7 + 4 pui'Ui 


(2.76) 


The resulting final form of the filtered LES energy equation is: 


A (£,) + £- fe£, + f, i P) 


d 

dx 3 


(lli 


T ij + U ‘ T ij ) 


dx 


(?, + «r) 


. i ( 


2.4 Turbulence Modeling 

Both the RANS and LES sets of equations derived in sections 2.2 and 2.3 require a 
turbulence model to close the momentum and energy equations. In the RANS 
approach, all unsteady turbulent motion is replaced by a turbulence model. The 
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resulting LES equations are very similar in appearance to the RANS equations, and 
also require a model to close the momentum and energy equations. The difference 
for the LES equations, however, is that the terms replaced hv a model are only the 
turbulent terms that are too small to be resolved using the filtered LES equations. 
As a result, the large scale turbulent motion is directly calculated, and the effects of 
the smallest scale turbulence are accounted for using a subgrid turbulence model. 

I he turbulence model employed here to close the RANS equations is the 
( ■ebeci-Smit.h algebraic turbulence model [18,19]. Since the RANS equations are 
only used in this hybrid method to calculate wall boundary layer regions with no 
adverse pressure gradients, the selection of a relatively simple algebraic model such 
as the Cebeci-Smith formulation is appropriate. The wall function technique of Ota 
and Goldberg [73] is used in conjunction with the Cebeci-Smith model to enable use 
of a computational grid with the first point off solid boundaries placed in the 
logarithmic layer. This wall function approach is based upon the compressible law 
of the wall formulation of White and Christoph [110, 111]. The filtered LES 
equations are closed using the Smagorinsky subgrid model [94]. 

Implementation of the wall function technique is critical to the development of this 
hybrid approach in order to enable use of a single computational grid extending 
continuously from the RANS regions to the LES regions. If a wall function 
approach were not used, grids for the RANS regions would have to be packed very 
tightly to the wall and use significant grid stretching, while a separate grid which 
minimizes grid stretching would need to be constructed for the LES region. Use of 
such non-continuous grids for the RANS and LES regions would require an 
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interpolation scheme that would likely introduce undesirable errors into the 
combined hybrid method. 

Use of the Cebeci-Smith model to close the RANS equations and the Smagorinsky 
subgrid model to close the LES equations is desirable in terms of code 
implementation. While the function of the Cebeci-Smith model to replace all ol the 
t urbulent stresses with a model is quite different from that of the Smagoi insky 
subgrid model, which only replaces the small subgrid turbulent stresses, both are 
eddv viscosity models and are derived at least in part, from mixing-length theory. 
The similar formulation of these two models enables the RANS equations and LES 
equations to be solved with a single solution scheme and computational giid. as 
mentioned previously. For a compressible nozzle or mixing layer flow, such as that 
depicted in figure 1.1. the change from RANS regions to LES region occurs at the 
vertical plane passing through the trailing edge of the splitter separating the wall 
bounded flows. 

In the following sections, details of the Cebeci-Smith turbulence model and the 
Smagorinsky subgrid model are provided. 

2.4.1 RANS Turbulence Model 

The unclosed terms from the RANS momentum and energy equations are the 
Reynolds stress shown in equation (2.28) and the turbulent heat flux, which is given 
in equation (2.40). As mentioned previously in section 1.1, the Boussinesq 
approximation is used to relate the turbulent Reynolds stress to the mean rate of 
strain tensor through a turbulent (or eddy) viscosity. This is directly analogous to 
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equation (2.9) which relates the viscous stress to the mean rate of strain through 
the molecular (or laminar) viscosity. The turbulent analogy to equation (2.9) is: 


T Ij = — /« 


Similarly, the turbulent heat 
turbulent conductivity. k T : 


= /' 




flux is related to the temperature gradient 


(2.78) 


through a 


q] = pu'p," 
c)T 


-k T 


(2.79) 


d.r :i 


I he turbulent Prandtl number, Pr ! is used to relate the turbulent viscosity to the 
turbulent conductivity: 


p r T 


l‘ T ('r 

k T 


(2.80) 


The turbulent Prandtl number is taken to be a constant here and equal to 0.9. 
Using equation (2.80) and assuming a constant turbulent Prandtl number enables 
the turbulent heat flux to be expressed as a function of the turbulent viscosity that 
is used to calculate the Reynolds stress. The turbulent heat flux becomes: 


% 


c P p T dr 


(2.81) 


Pr T dxj 

The Cebeci-Smith model, which was chosen here to close the RANS equations in the 
wall boundary layer regions, treats the wall boundary layer as having inner and 
outer regions where the turbulent viscosity is defined as: 
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y < fjm 


(2.82) 


\l- l luter- y > y ,n 

In equation (2.82). y,„ is defined as the smallest value of y (the distance away from a 
wall) at which The expressions for the inner and outer layer 

turbulent viscosities are as follows: 

Inner Layer: 



with the mixing length ( mtx is given by: 


(2.83) 





(2.84) 


Outer Layer: 


= op6 v U'F Uti (2.80) 

In equation (2.85). the quantity S v is the velocity thickness, u e is the boundary layer 
edge velocity, and F k ieb is the Klebanoff intermittency function. The velocity 
thickness is defined as: 


S v = 



( 2 . 86 ) 


This velocity thickness is identical to the displacement thickness for incompressible 


flows. 
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In references [52] and [53]. Klebanoff presents an (’xpression for ihe intermittency of 
turbulence near the edge of a boundary, which has a functional form involving the 
complimentary error function. This original intermittency function is usually 
approximated by the following formula, as indicated by Oebeci [18]: 


Huh = 


'+■«(!)' 


;i — 1 


(2.8 1 


The closure coefficients appearing in equations (2.84) and (2.85) are 


h = 0.40 


o = 0.0168 


.4+ = 26 


Wall Function Implementation: 

The Oebeci-Smith model is usually integrated down to the wall, using a 
computational grid with the first point off of the wall placed well within the laminar 
sublayer, corresponding to y + < 5. For the hybrid method developed in this work, 
the objective is to place the first point off of the wall in the logarithmic layer to 
enable the use of computational grids that are not packed as tightly to the wall. 
Removing the tight spacing requirement will enable a continuous grid into the LES 
region. In addition, because the allowable time step of the computations is 
proport ional to the size of the smallest grid cell, a less tightly packed grid enables a 
larger time step for the solution scheme. The wall function technique of Ota and 
Goldberg [73] is one of the more simple and effective methods currently in use, and 
it is the technique used in this work. 


Wall functions have been implemented most frequently in conjunction with 
two-equation k-e models. The benefits of implementing a wall function for use with 
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a k-< model are the same as that for the ( 'ebeci-Smith model used in this work 
including reducing grid requirements, and increasing the permissible time step of 
the computations. In addition, the need for the near-wall damping terms associated 
with low-Reynolds number k-c models is removed. These near-wall damping terms 
frequently result in numerical stiffness of the solution procedure. Nichols [i2] 
implemented the White-Christoph law of the wall with a k-c model for application 
to time dependent aerodynamic flows. Maui [64] implemented the Ota-Coldberg 
formulation in the WIND code for use with turbulence models ranging from 
algebraic to two-equation formulations. 

The use of a wall-function approach is st rictly only valid in flow regions absent of 
adverse pressure gradients and separations, due to the assumption that the law of 
the wall holds. However, Avva et al. [2] have shown results for separated flow’s in 
which wall function methods perform no worse than methods integrating to the 
wall. The intention of the wall function implementation in this w’ork is to only apply 
the method to attached wall boundary layers wdiere the law of the wall is valid. The 
wall function approach is compared to the standard procedure of integiating to the 
wall for two supersonic boundary layers in chapter 4. 

The Ota-Goldberg wall function employs the White-Christoph [110.111] 
compressible law of the wall: 



where the compressibility parameter 7 is given by: 
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(2.89) 


' 2CpT w 

In equation (2.88) uj is the value of u + at the first point off of the wall, y+ is the 
value of ir at the first point off of the wall, and y,} = 0.1287. In equation (2.89), 
the parameter r is the recovery factor, which is typically taken to be TV 3 for 
turbulent boundary layers, and T w is the wall temperature. An iteration procedure 
is used with equations (2.88) and (2.89) to solve for ut. from which the shear 
velocity a T can be obtained: 


= — (2.90) 

Finally. the shear velocity is used to compute the wall shear stress through: 

~r = f»<l (2.91) 

The wall shear stress calculated in equation (2.91) is then used in the solution 
scheme for the momentum and energy equations in the RANS regions. 

2.4.2 LES Subgrid Scale Model 

The terms that must be closed for LES equations are the subgrid-scale stress given 
by equation (2.28) and the subgrid scale heat flux, shown in equation (2.69). The 
earliest subgrid scale model for LES computations was developed by 
Smagorinsky [94]. Despite significant efforts to develop more sophisticated subgrid 
scale models, the Smagorinsky formulation is still widely used, and is itself the 
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foundation upon which some of the more sophisticated models are derived. 1 he 
form of the model is very similar to the Cebeci-, Smith model used for the R A XS 
equations, in that a gradient-diffusion mixing-length approach is used. 

T he Smagorinsky expression for the subgrid scale st ress is: 


^ = pfoiUj - 


= 2(Cs±) 2 p\/n {$ij 

The parameter tt is defined: 


2$kk-fiij) ~ §OA J /)ii S, 


(2.92) 


The parameter A is the filter width arid as a result, it is also used as the length 
scale that is characteristic of the subgrid turbulence. For use with a computational 
method, A is usually taken to be the grid spacing. In three dimensions, with a 
computational grid having unequal spacing in the three directions, this subgrid 
length scale is usually taken to be: 


A = ( A.rA;/A;-)' (2.94) 

For computational grids with substantially different spacing in the three directions, 
an alternative form (see Ragab [82]) is 


A = 


(A.r) 2 + (Ay) 2 + (Ar) 2 l ^ 
3 


(2.95) 
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In t lif' three dimensional simulations of a compressible mixing layer discussed later 
in this work, the two forms of the subgrid length scale presented in equations (2.9-1) 
and (2.95) are investigated. The constants C. s - and Cj have been found to be highly 
dependent on the flow under investigation. Rogallo and Moin [85] suggest a range 
for ( s in the range 0.10 < f .s 5: 0.24. Both of the limits on (’$ given b\' Rogallo 
and Moin are investigated for the mixing layer in this work. The constant ( ’/ is 
usually equal to 0.01. but several authors, including Ragab [82] and Choi et al. [20] 
mention that the contribution of the term involving Cj may not be important and 
may be neglected. This approach is taken in this work, and as a result, the original 
expression for the Smagorinsky subgrid scale stress in equation (2.92) may be 
rewritten as follows: 


, ~ d u, . 

; -o 3 /< dxj d,j 


where the subgrid scale turbulent viscosity is given by: 


(2.96) 


f*‘ 8 ‘ = p( ( 2.97 ) 

Note the similar form of equation (2.97) to the expression for the Cebeci-Smith 
inner region turbulent viscosity in equation (2.83). While the mixing length defined 
for the Cebeci-Smith model is used to characterize all of the turbulent motion, the 
length scale defined here for the Smagorinsky model only' characterizes the 
subgrid-scale motion. 

Finally, the subgrid scale heat flux is modeled analogously to that done for the 
turbulent heat flux of the BAN’S equations: 
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where k' 9 * is related to through the turbulent Prandtl number. As in the BANS 
regions, the turbulent Prandtl number is assumed to be constant in the LES regions 
and equal to 0.9. 1 he subgrid scale heat flux becomes 

< 9 ‘ _ ( _Z£Z 0L ( 2 . 99 ) 

lj Pr T ()s, 
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CHAPTER 3 


SOLUTION PROCEDURE 

The procedure used to solve the equations developed in chapter 2 is formulated in 
this chapter. The HANS and LES equation sets are recast in vector form, which 
corresponds directlv to the form of the equations that are solved by the 
computational method. The numerical scheme used to solve these equat ion sets is 
first illustrated using a one-dimensional model problem. The extension of the 
solution scheme to the RANS and LES equations in three dimensions on stretched, 
non-rectangular computational grids is presented next. The use of a generalized 
coordinate transformation, time-step selection, and artificial dissipation is discussed. 

3.1 Governing Equations 

Before discussing the numerical method used in this hybrid method, the RANS and 
and LES equations are expressed in vector notation, which combines the continuity, 
three momentum equations (corresponding to each of the cartesian coordinates), 
and the energy equation into a compact form that is actually used by the solution 
scheme. For reference, the RANS and LES equations derived using tensor notation 
in chapter 2 are repeated here. The equations written here differ from those in 
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chapter 2 only in that all terms now appear on the left side of the equations, which 
will make the conversion to the veclor notation more straightforward. 


In the HANS regions, the continuity, momentum, and energy equations are: 


Op 

dt 


+ 


o 

d.v t 


(3-1) 


d 0 

(p u i ) + 


dt 


Ox 


op Ot, 

(pu.llj) + 7-j -X— 

OX; OX; 


dxj 


(3.2) 


o_ 

dt 



~ ( ujE , + iijP) 


Ox 


0 , „ _ . T\ 0 

(">'o + «<”,) + 


()Xj 


(<Ij + <ij) = 0 


(3.3) 


Likewise in the LES region, the continuity, momentum, and energy equations are: 


c)p 0 
Ot Ox, 


= 0 


d_ 

dt 





. . dP 


dx 3 


dr-f 
- >J =0 


dx 3 


(3.5) 


d_ 

dt 



— (ujE' + ujP) 


0 

Ox, 


d 


+ UiT-f) + — (qj + q?') = 0 


Both the RANS equations (3.1 - 3.3) and the LES equations (3.4 - 3.6) can be 
expressed in Cartesian coordinates (x,v,z) using the compact vector form: 
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( 3 . 7 ) 


OG <)E OF OG 

— -| (- 1 

i:)t Ox Oy Oz 


For the RANS equations, the vectors Q, E. F, and G are: 


Q 


p 

pit 

pr 

prr 

Ei 


E 


pit 

pit 2 + P-t 
J yitv 


T xll ~ T 

xy 


F = 


puw — t x . 

a [E t + P] - u (r„ + rJJ - v [f jy + rjJ - w (r« + rJJ + y, + yj J 

pv 

pilV - T rll ■ 
pi 2 + P-T , 
pv w -T vz - T ‘. 

v [E t + P\ - u [t I y + rJJ - f (f M + r T J - te (r a; + r r J + <)„ + yj 


_T 


( 3 . 10 ) 


G = 


pic 

pine — t i: — rj. 


- T y*- T yz 

pic 2 + P — r_. — tJ, 

w [E t + P}- u (f x . + rJJ - r (r a2 + r r J - re (r„ + r 7 ! ) + <?_- + yj 


( 3 . 11 ) 
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In section 2.4.1. the form of the modeled t urbulent stress in equation (2.78) was 
shown to be the product of a turbulent viscosity and the mean rate of strain term. 

I his was directly analogous to the form of the viscous stress tensor shown in 
equation (2.9) which was defined as the product of a laminar viscosity and the same 
mean rate of strain term. This similarity allows the viscous and turbulent stresses in 
the flux vectors E, F. and G to be written together as the product of an effective 
viscosity. // T ft 1 and the mean rate of strain term. These combined stresses are as 
follows: 


- , 7 2 / T' i jj" fff The 

+ y = !<*.+ K r ) ( 2 $: - d " 


r.. + t; 


ij Ox 0 

Tx i > dw du c)v\ 


3 ( p + p )(^2— ----j 
T T r ( dil- dv 

'«» + T *y = + T y, = (/< + d ) [~r + 


V dy dx 
dw du \ 

«7 + fc) 

dr dw 


(- 1 - 12 ) 


T r, + Tj t = T ZI + t! x = (// + //' 

T - + T ‘‘ = f " + T “ = ^ + ,,T) { & + ly 

The same reasoning used to combine the viscous and turbulent stresses may be 

employed to combine the laminar and turbulent heat flux terms using an effective 
thermal conductivity k + k T . The combined heat flux terms in the vectors E. F. and 
G are 
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5, + 7 


5, 

q. + q r . = — (A* + />' f ) — 

For the LES equations, the vectors Q. E. F. and G are: 


()t _ 

( PpP 

+ 

Crp T 

\° T 

().r 

V Pr 

Pr T 

) ().r 

c)T _ 

(CpP 

1 

('pp T 

y± 

dy 

\ Pr 


Pr 1 , 

) dy 

<)t _ 

( Cph 

1 

CpP t 

ylL 

~d~z ~ 

\ Pr 

"T 

Pr T 

) dz 


Q 


p 

T>" 

pv 

pic 

E, 


(3.14) 



u [E t + P 


pit 

pu 2 + P - T xx - r;‘f 

pUV T xy 

pllW Tjr Z T xz 

(r. + C 9 ' ) - v {t xy + r;f ) - tc (T„ + t£‘) + <h + < 9t 

(3.15) 


pv 


t — r^- 

' xy xy 


F = 


/9(/r 

rtf? 2 4- P — T — T 39 * 

i n ' J yy 


pv* - T »-’ - r vf 


[£, + T 5 ] - u (r„ + T xy S ) ~ ~ v ( T yy + r ,T) “ l? ’ ( r »* + T w*) + ^ 


(3.16) 
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pCr 

pair — T r . — t r T 

G = prii' - t - r;f 

pi' 2 + P — t.. — r;T 

«’ [& + p] - II (t,, + r;f) - f (r y . + r’f) - tr (r,_. + pf) + <y ; + £»* 

(3.17) 

As was done for the RANS equations, the viscous and subgrid scale stresses 


appearing in the llux vectors E. F, and G are written together for the LES equations 
as the product of an effective viscosity p -f pr 9 ' and the resolved rate of strain term: 


. .... . . da dv (hi \ 


dv da die 

n, + c - +/'-> ^ 

.-> , , ( j)ir Oil dv \ 

t.,+t ;r = 3 , p + /r -..,( 2 -----) 


T xy T T~ g " — 7" yj, -|- 

1 jrj/ y j - ' yx 


(P + V 


du dv ' 
dy + c)x 


_ — , , / dir \ 

g- + r ,f = + 7-;, = + p ‘*’ ) ( 777 + 

( Or dw \ 

= n, + r; r = ( / , + ^-)(^ + ^) 


G-- + - - 


The laminar and turbulent heat fluxes are also expressed in a combined format 
using an effective thermal conductivity k -f k £gs \ 


<L + <1 


X 


% + <T 
+ qT 


~ {k + kS9!) lk = 

dT 

- (* + *•••) ^ = 


Cp/u Cpp’ 9 ‘ \ dT 

Pr Pr S9 -‘ / c/x 

C P p | C>f 

Pr Pr ,9S / 

Cp/z f '/•//-*' \ df 

Pr + p r *»* J dz 


(3.19) 
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3.2 Numerical Method 


The solution algorithm used in this work is the predictor-corrector scheme due to 
Gottlieb and Turkel [37]. This method was developed in the same philosophy as the 
original MacCormack scheme [61]. While the accuracy of the MacCormack scheme 
is second order in time and second order in space, the Gottlieb-1 urkel scheme is 
second order accurate in time and effectively fourth order accurate in space, and as 
a result , is often referred to as the "MacCormack 2—4 scheme. 

The Gottlieb- Turkel scheme has been applied to several time dependent flow 
problems because of its robustness, accuracy, and relative ease of implementation in 
methods for solving the Navier-Stokes equations. Snyder and Scott [95] 
demonstrated that the Gottlieb-Turkel was more accurate than other similar 
schemes for benchmark acoustic problems. Bayliss, et al. [8.9] successfully applied 
the scheme to boundary layer calculations. Ragab and Sheen [81,82] used the 
Gottlieb-Turkel scheme to perform LES calculations of compressible mixing layers in 
which the computational domain contained only the mixing region. 

The Gottlieb-Turkel scheme is illustrated using the following one-dimensional 
equation, which is a simplified model equation of the vector equation (3.7): 


The predictor step is: 


dt dx 


(3.20) 
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( t*i — ( h ~ ( — ' ft + 8/ , + 1 — J ,;+2 ) 


(3.21) 


The corrector stej) is: 


</” +i = ; 


At 




(3.22) 


The predictor-corrector scheme advances a solution in time from the time level (») 
to (■// + 1). 1 he time step. At , is related to the grid spacing, A.r. and the 
propagation speed, .4. through the CFL (Courant-Friedrichs-Lewy) number: 


At = CFL ~ (3.23) 

1 he Gottlieb- 1 Turkel scheme is cited by several authors, such as Hudson and 
Long [43], as providing second order accuracy in time and fourth order accuracy in 
space. Bayliss et al. [9] indicated that the scheme has fourth order accuracy only if 
At is of the order ( Ax) 2 , and Nelson [71] showed that the spatial accuracy of the 
scheme is only third order for CFL numbers approaching 1. A detailed analysis of 
the truncation error for the Gottlieb- Turkel scheme is provided in appendix A. This 
analysis shows that the leading truncation error terms resulting from discretization 
of the model equation 3.20 is obtained from the following equation: 

()q ( Qf At 2 cFq Ax 3 r) 4 f Ax 4 c)\f 

dt + Ox ~ 6 OP L 18 dx* + 30 IhC (3 ' 24) 

The left side of this equation is the form of the model problem in equation (3.20) 

while the right side is the truncation error. The first term on the right side of 
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equation (A. 13) indicates that the scheme provides second order accuracy in time. 
The second truncation error term indicates that the scheme is only third order 
accurate in space for CFL numbers approaching I. as previously shown by 
Nelson [71]. However, for most problems in which the Gottlieb-Turkel scheme is 
used, the maximum CFL number is usually set to a value of O.o or less. In addition, 
for computational grids that employ non-uniform stretched grids, the limiting time 
step is inversely proportional to the smallest grid spacing and the resulting effective 
CFL number will be much smaller in regions where the grid spacing is larger. This 
is demonst rated by considering equation (3.23) for the case of vaiiable giid spacing. 
A.r, but for a constant time step \t. For such regions of the computational domain, 
the second truncation error term on the right side of equation (3.24) will be 
insignificant, and then the next truncation error term indicates that the scheme 
provides fourth order spatial accuracy. In conclusion, the Gottlieb-Turkel scheme is 
strictly second order accurate in time and third order accurate in space, but in the 
case of computational grids with significant stretching, the spatial accuracy is 
effectively fourth order for most of the computational domain. 

3.3 Extension of the Gottlieb-Turkel Scheme to 
Generalized Coordinates 

The Gottlieb-Turkel predictor-corrector scheme is extended to three dimensions and 
generalized curvilinear coordinates in this section. The use of generalized 
coordinates, necessary for the solution of the RANS and LES equations on 
computational grids with non-uniform stretched spacing and non-iectangulai gtid 
cells, is presented in section 3.3.1. The time step calculation procedure is discussed 
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in section 3.3.2 and the use of numerical dissipation for stability, particularly in the 
HANS regions, is discussed in section 3.3.3. 

3.3.1 Generalized Coordinates 

The three dimensional vector forms of the HANS and LES equations, both given by- 
equation (3.7). are an extension of the one dimensional problem given in equation 
(3.20). The Gottlieb- 1 urkel scheme can be written for three dimensional problems 
as follows; 


Q +1 = C : CyC T P~P y P J .Q n (3.25) 

where P r . P, r and P- are the one-dimensional predictor operators and C r . C y . and 
('• are the one-dimensional corrector operators that correspond to the flux vectors 
E. F. and G respectively. 

The numerical scheme is further extended to generalized coordinates through the 
transformation of the RANS and LES equations from cartesian physical space 
[x.y.z) to computational space (£. //. 0- Hixon et. al. [41] found that the chain rule 
formulation was more accurate than other formulations using generalized curvilinear 
coordinates. The chain rule formulation is used in this work, in which the flux 
vector derivatives are expressed as: 


0E 0E OE . dE 
lh- ~ ix ~d^ + T1r lhi + L,I 'd^ 

OF dE dE dE 

-t'W + + ''’di : 

dG r dG dG . dG 
dz d$ +7h ~d^ + (,: l^ 


(3.26) 
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The terms »/,. C* C»- 'h- all<1 C are the metrics of the transformation 

from physical space to computational space. The derivation of these metric terms 
are detailed in appendix B. They are computed using a fourth order finite diffeieiue 
method to be consistent with the (lottlieb-Turkel scheme. 

The stress and heat flux terms appearing within the flux vector expressions for the 
RANS and LES equations were shown to involve the derivatives of the velocity- 
components and temperature with respect to the cartesian coordinates ,c. y. and ~ 
in equations (3.12. 3.13) arid (3.18. 3.19). These velocity and temperature 
derivatives are also computed using the chain rule form. For example, the HANS 
derivative of the mean velocity u with respect to x is: 


du , du da . du 

7- = + >h: 7— + <,,7T7 

dx di ay ()(, 


(3.27 


Similarly, the LES derivative of the resolved velocity u with respect to x is: 


dii 

Ox 


d u d u . d ii 

+ "-a7 + l ”ac 


(3.28) 


To compute the derivatives of the velocity and temperature terms with respect to 
the computational coordinates £. y. (,'■ Bayliss et. al. [9] have shown that a specific 
procedure is required retain the overall accuracy of the Oottlieb Tuikel scheme. Foi 
the predictor step, which uses forward differencing for the flux terms as shown in 
equation (3.21), the velocity and temperature derivatives taken in the same 
computational direction as the flux vector derivative are obtained with a two point 
backward difference operation, while the derivatives taken in the other two 
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directions are obtained with a central difference operation. For example, the u 
velocity derivatives taken to compute the flux vector terms ^ . and in 
equation (3.26) are: 


a, - 


c)ii 

dit i A 

Du 1 . 


( :i.29 ) 


<K 


- \ Uk + 1 — "k-l , 


A similar procedure is used for the corrector step, for which the flux terms are 
backward differenced, as shown in equation (3.22). The velocity and temperature 
derivatives taken in the same computational direction as the flux are calculated 
using a two- point forward differenced difference expression, while the derivatives in 
the other two directions are calculated with a central difference operation. 


3.3.2 Time Step Calculation 


A procedure for calculating a time step that will allow for stable calculations with 
predictor-corrector schemes was proposed by MacCormack in reference [62], The 
procedure involves searching the entire computational domain for the minimum 
value of: 


A/ 


| u | | v | | w | I l 1 I \ 

A.r A;y + Ac + V (Ax ) 2 ' (A //) 2 + (Ac) 2 / 


(3.30) 


The minimum time step obtained from equation (3.30) was modified through the 
use of the CFL number, which was always less than one for the Gottlieb-Turkel 
scheme, so that the actual time step used in the computations is given by: 
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A = CFLM 


(3.31) 


For the hybrid calculat ions of the compressible mixing layer, the calculat ed time 
step is allowed to vary during the initial transient part of the calculation, when the 
mixing layer is allowed to develop. For the rest of the calculation, however, the time 
step is fixed to a constant value corresponding to the minimum value observed 
during the initial transient development . A constant time step is necessary for the 
LES region to calculate turbulent statistics from flowfield data stored at constant 
intervals in time. 

3.3.3 Numerical Dissipation 

Numerical dissipation is usually required in Euler and RANS computations to 
remove numerical oscillations that are undesirable in terms of solution accuracy and 
code stability. These oscillations are typically at large wave numbers caused by 
nonlinearities in the solution process. Jameson et al. [46] and Pulliam [80] have 
developed numerical dissipation schemes that have been used extensively for Euler 
and Navier-Stokes calculations of aerodynamic flows. These schemes operate by 
adding dissipative terms to the equations of motion, and require careful use to keep 
the levels of numerical dissipation as small as possible while retaining stability of 
the solution. 

For LES computations, the filtering process is used to develop a set of equations 
that, only resolve the large scale unsteadiness associated with smaller wave numbers. 
The larger w r ave number unsteadiness, including both the subgrid scale turbulence 
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and the unphysical numerical oscillations, are handled through the use of a subgrid 
scale model. 1’he eddy viscosity formulation of the Smagorinskv model used in this 
work effectively adds a viscous stress to the momentum and energy equations which 
serves to both replace the unresolved subgricl scale stresses and to damp unphysical 
oscillations away from the flow regions of interest. 

Some authors such as Kennedy and Carpenter [48.49] and Hixon [40] have 
developed explicit filters for damping unphysical large wave number oscillat ions in 
time dependent flow calculations. These filters are not directly associated with t he 
filtering process, described previously, used to derive the LES equations, but instead 
are another class of numerical dissipation schemes. Hixon [40] has shown that the 
explicit filters of Kennedy and Carpenter may be used in conjunction with explicit 
solvers such as the Gottlieb-Turkel scheme to damp unresolved oscillations, and as a 
result, this method is employed here. 

The Kennedy and Carpenter filters have the form: 


Q = (1 + o n D) Q (3.32) 

where D is a symmetric matrix filter function of order 2 n(n = 1,2,3....) and the 
coefficient <>,, is given by: 


-(!)" 


(3.33) 


Kennedy and Carpenter developed the filter o n D to be of a form that retains larger 
wave number components with increasing n. and to always be dissipative in effect. 
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The matrix function D and the coefficient o„ are derived for filters of orders 2 
through 20. corresponding to n = 1. 10 in reference [40]. In this work, an eighth 
order filter (v = 4) is used, which has the form in one dimension: 

( — Qi —4 + $();-:< - ‘ddc /,-2 + - 70 qi + 56c/, +, - 2 Sc/, + 2 + 8r/, +:i - </ l+ 4 ) 

°» D = 256 

(3.34) 

Skewed stencils of similar form to equation (3.34) are used for points at and near 
boundaries. The filter function is applied everywhere in the HANS region, but is 
only applied in the LES region when a transient in pressure indicates that the 
solution scheme will numerically become unstable. As a result , the procedure used 
in the LES region is to modify the coefficient o n to be a function of the local 
pressure, such that when the pressure at a point m the flow drops lowei than 10 
percent, of a reference pressure characteristic of the flow, o„. is modified to be 


_ -(1)” ( ■_ IQ/VA 4 
Qn " (2) _2n V P ) 


(3.35) 
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CHAPTER 4 


RANS CALCULATIONS OF WALL BOUNDARY LAYERS 

In this chapter, a series of flat plate boundary layers are calculated to investigate 
the RANS solver of the hybrid method. Since no free shear layer mixing regions 
occur in these cases, the LES part of the method is not used in these calculations. 
All of these boundary layer calculations were obtained for a. flowfield situation 
depicted by figure 4.1 where a uniform inflow meets a smooth flat plate and a 
boundary layer develops over the plate surface. 

The first case considered here is an incompressible laminar boundary layer that is 
calculated using the RANS solver without any turbulence model. The second case is 
an incompressible turbulent boundary layer that is calculated with the 
Cebeci-Smith turbulence model integrated to the wall, without the use the wall 
function approach. This will be referred to hereafter as the “wall-integration 
approach. The third and fourth cases investigate two supersonic boundary layers 
that have the same flow' conditions as the two isolated flow's that form the mixing 
layer w-hich will be investigated in the next two chapters using the combined 
RANS-LES procedure. These last two boundary layer calculations are run with 
both the wall-integration and wall-function approaches. With the formulation of the 
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Uniform 

Inflow 



Flat Plate 


Figure 4.1: Schematic of flat plate boundary layer 


RANS solution procedure to march in time using the Gottlieb- Turkel scheme, the 
flowfields for all of these cases are initially set to the freest ream conditions 
everywhere in the computational domain, and the calculations progress until a 
steady state solution is obtained. The residual error and flowfield quantities such 
the skin friction coefficient were used to monitor convergence to steady state. 
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Figure 4.2: Computational grid for the incompressible laminar flat plate case 

4.1 Incompressible Laminar Boundary Layer 

The flow used for the laminar boundary layer case was constructed such that the 
Reynolds number was 10000 at the end of the plate, using the freestream density, 
velocity, and viscosity. The vertical dimension of the flow domain was one half that 
of the axial dimension. Because the flow solver was developed using the 
compressible form of the Navier-Stokes equations, the freestream flow was set to 
correspond to Mach 0.2. Use of a smaller freestream Mach number would 
substantially slow down the convergence characteristics of the numerical scheme. 

Figure 4.2 shows the computational grid used for the calculations with 51 axial 
points bv 51 vertical points. Calculations were initially obtained using this two 
dimensional grid and a two-dimensional version of the flow solver. Additional 
calculations were obtained with the three dimensional solver to fui ther validate the 
implementation of the computational method. For the three-dimensional 
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calculations, 11 points were used in the spanvvise direction to accommodate the 
three-d imensional flow solver, hut no flow development occurs in this direction. The 
grid was packed to the leading edge in the axial direction such that the grid spacing 
from the first to second point corresponds to a Reynolds number of 10, again using 
freest ream flow properties. In the vertical direction, the grid was packed to the plate 
surface such that the grid spacing from the first, to second point corresponds to a 
Reynolds number of 2. A hyperbolic tangent stretching function was used to stretch 
the grid both in the axial and vertical directions. For the three-dimensional 
calculations, the grid points were equally spaced in the spanvvise direction. 

The boundary condition for the flat plate is set to be a no-slip, adiabatic surface 
such that all of the velocity components are set to zero, and the vert ical 
temperature gradient is zero. At the inflow, the total pressure and total 
temperature are specified to correspond to Mach 0.2 flow at sea level atmospheric 
conditions. The static pressure is set at the outflow to also correspond to Mach 0.2 
flow at sea level atmospheric conditions. The boundary along the plane at the 
highest vertical point is modeled as a slip-wall surface, which is also known as a 
symmetry surface. A boundary condition in which all flow quantities were 
extrapolated from the interior was also investigated for the top boundary, and 
produced identical results for that obtained with the slip surface. For the three 
dimensional calculations, the boundaries along the two extreme spanwise planes 
were also set as slip surfaces. Because the solutions obtained with the two and three 
dimensional calculations were found to be nearly identical, only one set of computed 
results are discussed in the the rest of this section. 
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The exact solution for the incompressible flat plate boundary layer with constant 
fluid properties and no pressure gradient was obtained by Blasius [11]. and as a 
result, is known as the Blasius solution. The Blasius solution enables velocity 
profiles at different axial positions, ,r to be reduced to a single similarity profile, 
through the use of a similarity variable. ?/. defined as: 


V = 




The exact solution of Blasius enables the skin friction coefficient, boundary layer 
thickness, displacement thickness, and momentum thickness to be expressed as a 
funct ion of axial position. 

Skin friction coefficient: 


0.664 

Cf ~ \7C 

Boundary layer thickness: 


6 _ 5 

X \fR<~v 

Displacement thickness: 


5 * 1.7208 

7 " 

Momentum thickness: 


( 4 . 2 ) 


( 4 . 3 ) 


( 4 . 4 ) 
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figure 4.3: Skill friction for incompressible laminar boundary layer 


0 _ 0.664 

x “ v/7?v; 


(4.5) 


The calculated skin friction along the flat plate is compared to the expression from 
the Blasius solution given in equation (4.2) in figure 4.3. The overall agreement 
between the calculation and the expression from the Blasius solution is good. At the 
end of the plate, near Rt r ~ 10000, the calculation indicates that the skin friction 
levels off. This is a numerical effect of the outflow boundary condition which 
extrapolates the velocity components at the outflow plane from the interior, and 
indicates that the outflow plane should be downstream of the flow region of interest 
when such an outflow boundary condition is used. 
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Comparisons of the calculated velocity profiles are made to the exact solution ol 
Blasius in figure 4.4 at two axial positions corresponding to R< = 2000 and 

r - 5000. Again, the agreement between t he exact solution and the calculated 
results is close. Comparisons of the calculated boundary layer, displacement, and 
momentum thicknesses to those obtained from the Blasius solution are shown in 
figures 4.5, 4.6, and 4.< respectively. 

The boundary layer thickness along the flat plate obtained from the calculations was 
determined to be the distance from the wall where the local axial velocity became 
99 percent of the freestream velocity. The displacement and momentum thicknesses 
were obtained by numerical integration using equations (4.6) and (4.7) lespectiveh . 


<T 




(4.6) 


e = 




(4.7) 


These forms of the displacement and momentum thicknesses are valid for both 
laminar and turbulent incompressible flows. For compressible flows, the variation in 
density appears in the two expressions, which will be shown later in the discussion 
on the compressible boundary layer calculations. 


The agreement between the exact Blasius solution and the computed solutions 
obtained with the Gottlieb-Turkel scheme is very close for the boundary layer 
thickness and the integral thicknesses. At the end of the plate, A/I = T all of the 
three thicknesses flatten in slope. As was the case for the skin friction comparison in 
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Figure 4.4: Velocity profiles for incompressible laminar boundary layer 
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Figure 4.5: Boundary layer thickness for incompressible laminar boundary layer 



Figure 4.6: Displacement thickness for incompressible laminar boundary layer 
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I' igure 4.7: Momentum thickness for incompressible laminar boundary layer 


figure 1.3. this is due to the outflow boundary condition which extrapolates the 
velocity components from one point interior to the computational grid. 

4.2 Incompressible Turbulent Boundary Layer 

The second boundary layer case examined was an incompressible turbulent flow over 
a flat plate at Mach 0.2. A computational grid having 71 axial points by 71 vertical 
points was examined with the two-dimensional flow solver. The computational 
domain extended to a position corresponding to a Reynolds number of 2,000.000. 
The vertical dimension was equal to 40 percent of the axial dimension. The 
computational grid, shown in figure 4.8, was packed to the flat plate surface such 
that the average y + of the first point off the wall was set to 2.5. This average y + 
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that was employed to construct the computational grid was determined 1>v the 
following procedure: 

The definition of the wall normal coordinate t/ + is: 


+ V u r 
V = — 7 — 
Pw/Pw 

Expanding the expression for the shear velocity: 


(4.8) 


« T = 




(4.9) 


The expression for (y + ) in equation 4.8 can be rewritten in terms of the local skin 
friction coefficient and flow properties at the wall and in the freestreani: 


y+ = ( fKt^) ( 4 . 10 ) 

y V 2 //, x VV P-< Pn. / 

For this incompressible flat plate boundary layer case, the flow' properties are 
constant through the flow' such that p w = poo and y w = p 0 ©• Using a skin friction 
coefficient that is characteristic of the flown in this case C / ~ 0.003, and the 
freestream conditions, equation (4.10) can be used to calculate a characteristic or 
"average" y + for the grid spacing at the wall. In the axial direction, the 
computational grid was packed to the leading edge such that the initial spacing was 
A.r + = 75 using the same method used to calculate the average y + . A hyperbolic 
tangent stretching function w'as used to stretch the grid both in the axial and 
vertical directions, as was done for the laminar flat plate case. 
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Figure 4.8: Computational grid for t he incompressible turbulent flat plate case 


Thi s incompressible turbulent boundary layer flow was calculated using the same 
HANS approach as that for the laminar boundary layer, with identical boundary 
conditions, except that the Cebeci-Smith turbulence model was also employed. The 
wall function approach was not used here. With the grid packed tightly to the plate 
surface and the first point off the wall placed well within the laminar sublayer, the 
Cebeci-Smith turbulence model was integrated to the wall. 

Unlike the laminar boundary layer, an exact solution is not available for this 
turbulent flow. As a result , the computed results are compared to a benchmark 
experimental data set of Wieghardt and Tillman [112] for the wall skin friction 
coefficient and velocity profiles. Correlations for the boundary layer, displacement 
and momentum thicknesses obtained by Schlichting [90] are used for comparison to 
the computed results. 
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F igure 4.9: Skin friction for the incompressible turbulent boundary layer 


A comparison of the calculated skin friction with the data of \\ ieghaidt and 
Tillman is provided in figure 4.9. The calculated velocity profile at a location along 
the plate corresponding to a Reynolds number based on axial position of 1,050,000 
is compared to experimental measurements in figure 4.10. Figure 4.11 pi ov ides a 
comparison of the same results for the velocity profile using wall coordinates. 
Overall, the agreement between the calculated boundary layer results and the 
experimental data of Wieghardt and Tillman is very close. 

A comparison of the calculated boundary layer, displacement, and momentum 
thicknesses to correlations provided by Schlichting [90] are shown in figures 4.12, 
4.13, and 4.14 respectively. The correlations for these thicknesses are as follows 
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Figure 4. 10: Velocity profile for the incompressible turbulent boundary layer at Re x 
1.050.000 
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Figure 4.11: Velocity profile for the incompressible turbulent boundary layer at Re r 
1,050,000 using wall coordinates 
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Figure 4.12: Boundary layer thickness for incompressible turbulent boundary layer 



Figure 4.13: Displacement thickness for the incompressible turbulent boundary laye 
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Fi gure 4.14: Momentum thickness for the incompressible turbulent, boundary layer 


Boundary layer thickness: 


S -i 

- = 0.37 Re r s (4.11) 


Displacement thickness: 



(4.12) 


Momentum thickness: 


9 = 



(4.13) 
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The correlations provided in equations (4.11 - 4.14) are derived assuming the 
velocity profile for an incompressible turbulent flow over a smooth surface can be 
represented by a =-th power law profile. That is: 


The agreement between the correlations and the computed solutions obtained with 
the Gottlieb- Turkel scheme and Cebeci-Smith turbulence model is good. Again near 
the end of t he plate, the outflow boundary condition extrapolates the components of 
the velocities such that the three thicknesses flatten between the next to last giid 
point in the axial direction and the outflow boundary. This behavior does not 
adversely affect the rest of the boundary layer development upstream of the outflow 
boundary as indicated by figures 4.12 - 4.14. 

4.3 Compressible Turbulent Boundary Layers 

In this section, two compressible supersonic boundary layers are investigated to 
further validate the RANS approach. The two boundary layers selected were the 
two flows entering the mixing layer to be examined with the hybrid RANS/LES 
method. The particular mixing layer is one of the experiments of Goebel and 
Dutton [34-36] referred to as Case 2 in these references. Their experimental 
configuration is depicted by the schematic shown in figure 1.1. Upstream of the 
mixing layer, the flow on the top half of the splitter plate is for a Mach 1.91 stream 
and a total temperature of 578 I\. On the bottom half of the splitter plate, a 
boundary layer grows for a Mach 1.36 stream at a total temperature of 295 K. The 
Goebel-Dutton experiments provide details of the boundary layer entering the 


NASA/TM — 2001-210811 


81 



mixing region including measurements of the boundary layer, displacement, and 
momentum thicknesses. More details of the experimental configuration, flow 
conditions, and data collection, particularly regarding the mixing section, will be 
provided in chapter ■'3. 

For each of the supersonic boundary layers, both the wall-integration and 
wall-function approaches are investigated. As mentioned in chapter 2. the 
motivation for employing the wall function approach in conjunction with the 
Cebeci-Smit.h turbulence model was to enable a continuous computational grid from 
the wall bounded HANS regions into the mixing layer that is then calculated using 
the LES part of the hybrid approach. Comparison of the wall-function approach to 
the more standard wall-integration approach and the compressible boundary layer 
correlations determines the capability of the wall-function implementation to 
accurately provide the mean flow characteristics of the incoming boundary layer. 

In the previous section, calculated results were compared to benchmark 
experimental data for the incompressible turbulent boundary layer. Although the 
measurements of the boundary layer, displacement thickness, and momentum 
thicknesses at the trailing edge of the wall bounded regions in the Goebel-Dutton 
experiments make them one of the better documented mixing layer data sets, other 
details of the boundary layer developments are not provided. 

As a result, well established correlations for the skin friction coefficient, velocity 
profiles, boundary layer thickness, and integral thicknesses are used to investigate 
the details of the boundary layer development for the two particular supersonic 
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flows here. The calculated skin friction coefficient is compared to a combined 
correlation, that, is referred to here as the Mager-Schlichting correlation. In 
reference [90]. Schlichting provided the following expression for the skin friction 
coefficient for an incompressible flat plate with turbulent flow from the leading edge. 


Cj = .0592 /iV x * (4.15) 

In reference [63], Mager provides a correction to equation (4.15) for the skin friction 
coefficient in compressible turbulent flow: 


/ nnc. 

L / 



(4.16) 


To evaluate the other boundary layer characteristics such as velocity profiles, 
boundary layer thickness, and the integral thicknesses, the method developed by 
Tucker [105] is used here. Tucker provides a power-law relationship for the velocity- 
profile in turbulent compressible flow that is similar to the incompressible form 
given in equation 4.14. Tnlike the incompressible form which assumes a constant 
exponent of 7. the compressible form uses a variable exponent. 





where the exponent N is given by: 


(4.17) 


N - ( Re a m ) U 


(4.18) 
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I he Reynolds number lit (i m is based on the length of t he boundary layer 
development. The subscript am refers to the procedure of evaluating all fluid 
properties at the arithmetic mean of the wall temperature and the freest ream 
temperature. Tucker shows that a solution of the karman momentum equation 
results in: 


The parameter K is given by: 

A = 0.0131 ( ) (4.20) 

V Po(i 0 / 

where ft 0 , p 0 . and a 0 are the dynamic viscosity, density, and speed of sound, 
respectively, all evaluated at the freestream total temperature. 

Finally, the Mach number parameter, m 2 . is given by: 

2 M 2 

m 2 = — - (4.21) 

Tucker provides further integral relationships between the momentum thickness. 0, 
given by equation 4.19 and the boundary layer thickness, S. the displacement 
thickness, 5*. While the ratios 0/S and 8*/S are assumed to be nearly constant for 
incompressible flows, as indicated by equations 4.12 and 4.13, Tucker shows them to 
be a function of the freestream Mach number and the velocity profile parameter, N. 
Reference [105] provides extensive tables from which S. S *, and 0 may be obtained. 


(1 + m 2 ) 2 

U'(i + ffJ 


x ■ 


(4.19) 
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Following the validation of the HAN'S approach for the Mach 1.91 and Mach 1.36 
cases, the objectives will be to provide boundary layers to the mixing section in the 
hybrid RANS/LES calculations that match the mean flow properties measured in 
the Dutton-Goebel experiment. As a result, physical dimensions will be used to 
describe the computational grids and the calculated thicknesses 3, 3 , and 0 in the 
following two sections. 

4.3.1 Mach 1.91 Flat Plate Flow 

The computational grid for the wall-integration case is shown in figure 4.15 and the 
grid for the wall function case is shown in figure 4.16. Both grids extended 300 mm 
in the axial direction and 150 mvn in the vertical direction. In addition, both giids 
used 141 points in the axial direction arid 141 points in the vertical direction. The 
wall-integration grid was packed to the wall such that the first grid spacing was 
0.006 mm, corresponding to an average y + of 2.5, using equation (4.8). The wall 
function grid had the first point placed at 0.05 mm. corresponding to an average y + 
of approximately 20. This wall spacing was chosen for use with the wall-function 
grid so that the initial grid spacing at the wall was exactly l/10th of the splitter 
plate thickness in the experiment. This grid spacing could then be continued into 
the mixing region, with 10 points spaced equally in the vertical direction at the base 
of the splitter. In the axial direction, both the wall-integration and wall-function 
grids were packed to leading and trailing edges with spacings at the two ends set to 
0.10 mm, corresponding to l/5th of the splitter base thickness. 

The wall boundary condition for the calculations was set to be an adiabatic no-slip 
surface and the extreme vertical boundary was set as a slip surface, as was done for 
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Figure 4.15: Computational grid for the Mach 1.91 turbulent flat plate case using the 
vval 1-i ntegrat ion met hod 



Figure 4.16: Computational grid for the turbulent Mach 1.91 flat plate case using the 
Ota-Goldberg wall function 
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the incompressible turbulent boundary layer. For this supersonic boundary layer, 
however, the inflow condition was fixed and the outflow condition extrapolated 
quant ities. 

The skin friction coefficient along the flat plate obtained from the two computed 
solutions is compared to the Mager-Schlichting correlation m figuie 4.1 <• The 
wall-integration skin friction results were obtained by directly evaluating t he shear 
stress at the wall, since the grid was packed to the wall surface with the hist point 
off the wall placed in the laminar sublayer. For the wall-function case, however, the 
shear velocity that is solved with an iterative procedure within the Ota-Cioldberg 
wall-function was used to calculate the skin friction coefficient, along the flat plate. 
The agreement of the two solutions with the Mager-Schlichting correlation is good, 
with the wall-integration approach providing closer agreement. 

The calculated velocity profiles at a location along the plate corresponding to a 
Reynolds number based on axial position of 4.000,000 is compared to results of the 
Tucker method in figure 4.18. The profile obtained with the Tucker correlation was 
obtained by using equation (4.17) with the exponent N given in equation (4.18) 
evaluated to be 6.275 at. this axial position. In figure 4.19 the velocity profile at the 
same location expressed in w r all coordinates, and y + is compared to the 
White-Christoph compressible law of the wall given in equation (2.88). Figures 4.18 
and 4.19 indicate that both the wall-integration and wall-function calculations 
provide good agreement with analytical expressions developed to evaluate 
compressible turbulent boundary layer characteristics. 
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Figure 4.17: Skin friction for the Mach 1.91 turbulent boundary layer 



0 0.2 0.4 0.6 0.8 1 


u/u 

X 


Figure 4.18: Velocity profile for the Mach 1.91 turbulent boundary layer at Re 
4.000.000 
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Figure 4.19: Velocity profile for the Mach 1.91 turbulent boundary layer at R( r — 
4.000.000 using wall coordinates 


The calculated boundary layer, displacement, and momentum thicknesses are 
compared to results obtained using the Tucker method in figures 4.20, 4.21, and 
4.22 respectively. The boundary layer thickness along the turbulent compressible 
flat plate was determined in the same manner used for the incompressible tuibulent 
case, by finding the distance from the wall where the local axial velocity became 99 
percent of the freestream velocity. The displacement and momentum thicknesses for 
the compressible case differ from the incompressible case in that the density is 
involved in the integral expressions. These two integral quantities are defined for the 
compressible case in equations (4.22) and (4.23): 
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Figure 4.20: Boundary layer thickness for the Mach 1.91 turbulent boundary layer 
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Figure 4.21: Displacement thickness for the Mach 1.91 turbulent boundary layer 
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Figure 4.22: Momentum thickness for the Mach 1.91 turbulent boundary layer 
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The three measures of the boundary layer growth shown in figures 4.20 - 4.22 again 
indicate close agreement with results obtained using the Tucker analysis. The 
greatest discrepancy in these results is for the wall-function solution neai the 
leading edge of the plate. This is a result of the wall-function grid having 
significantly fewer points than the wall- integration grid near the wall to resolve the 
thin boundary layer at the leading edge. Further downstream, however, the 
wall-function approach provides similar agreement to the Tucker correlation. 
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4.3.2 Macli 1.36 Flat Plate Flow 


The computational grids used for the wall-integration and wall fund ion cases are 
shown in figures 4.23 and 4.21 respect ively. The construction of these grids was 
similar to that performed for the Mach 1.91 calculations, with the exception of the 
axial dimension being shorter for these Mach 1.36 calculations. This was done 
bt ’cause the Mach 1.36 boundary layer thickness m the Dutton-Ooebel experiment 
measured at the beginning of the mixing section was smaller than the Mach 1.91 
boundary layer, as will be discussed in the next chapter. The grids extended 200 
nun in the axial direction and 150 nun in the vertical direction. Both grids used 141 
points in the axial direction by 141 points in the vertical direction. The grid for the 
wall-integration case was packed to the wall such that the first grid spacing was 
0.006 nun. corresponding to an average y + of approximately 3.0. again using 
equation (4.8). The wall function grid had the first point placed at 0.05 nun. 
corresponding to an average y + of approximately 25. Both grids were packed in the 
axial direction to leading and trailing edges with spacings at the two ends set to 
0.10 nun. corresponding to 1 /5th of the splitter base thickness. The boundary 
conditions used for these Mach 1.36 cases were the same as those used for the Mach 
1.91 case, with the only difference being that a Mach 1.36 flow was fixed at the 
inflow here. 

The skin friction coefficient obtained from the two computed solutions is compared 
to the Mager-Schlichting correlation in figure 4.25. The agreement of the two 
solutions with the correlation is nearly as close as that for the Mach 1.91 case, with 
the directly calculated shear stress from the wall-integration approach providing 
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Figure 4.23: Computational grid for the turbulent Mach 1.36 flat plate case using the 
wall-integration method 



Figure 4.24: Computational grid for the turbulent Mach 1.36 flat plate case using the 
Ota-Goldberg wall function 
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Figure 4.25: Skin friction for the Mach 1.36 turbulent boundary layer 


closer agreement to the correlation than the wall-function approach. The 
calculated velocity profiles at a location along the plate corresponding to a Reynolds 
number based on axial position of 4.000,000 is compared to results of the Tucker 
method in figure 4.26. For the Mach 1.36 boundary layer at this location, the 
Tucker analysis indicates that the velocity profile exponent is 6.370. In figure 4.27 
the velocity profile expressed in wall coordinates, u + and y + is compared to the 
White-Christoph compressible law of the wall profile. In figures 4.26 and 4.27, both 
the wall-integration and wall-function calculations indicate acceptable agreement 
with the analytical expressions. 

Finally, the calculated boundary layer, displacement, and momentum thicknesses 
are compared to the Tucker method in figures 4.28, 4.29, and 4.30 respectively. The 
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Figure 4.26: Velocity profile for the Mach 1.36 turbulent boundary layer at Rt x 
4. 000, 000 



Figure 4.27: Velocity profile for the Mach 1.36 turbulent boundary layer at Re r 
4, 000. 000 using wall coordinates 
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Figure 4.28: Boundary layer thickness for the Mach 1.46 turbulent boundary layer 



agreement of the solutions with the Tucker correlation is good, although the 
agreement is not as close as for the Mach 1.91 case. The wall-function solution 
again indicates the largest discrepancy near the leading edge of the plate where the 
number of grid points available to resolve the thin boundary layer is low. The 
discrepancy is minimized further downstream, where the boundary layer becomes 
thicker. 
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Figure 4.29: Displacement thickness for the Mach 1.36 turbulent boundary layer 
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Figure 4.30: Momentum thickness for the Mach 1.36 turbulent boundary layer 
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CHAPTER 5 


TWO DIMENSIONAL MIXING LAYER CALCULATIONS 

Two dimensional calculations with the hybrid RANS-LES method were performed 
for a benchmark compressible mixing layer experiment and the results are described 
in this chapter. While true LES simulations require computations in three spatial 
directions, it is useful to compare two dimensional calculations to investigate effects 
of the RANS-LES interface region, axial grid resolution, and boundary conditions. 
The results of three dimensional calculations are presented in chapter 6 for the same 
experimental configuration. 

In this chapter, details of the experimental configuration and operating conditions 
are provided in section 5.1. Construction of the computational model is provided in 
5.2. Two-dimensional calculations investigating axial grid density effects are 
presented in section 5.3. Finally, calculations which investigate the splitter plate 
thickness and mixing section wall placement are presented in section 5.4. 

5.1 Experimental Configuration 

The flow' that is the focus of the hybrid method calculations is a two-stream, 
turbulent planar mixing layer that was examined in the experiments of Goebel and 
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Dut t.on [:M ■:{()]. A simplified schematic of their experimental configuration provided 
in figure o.l shows that two isolated air streams, in which boundary layers develop 
over a splitter plate surface, are brought together into a constant area mixing 
section. In all of their experiments, the higher speed primary stream occurred over 
1 1n' top surface of the splitter plate. 1 lie top stream enters the mixing sect, ion 
axially, while the bottom stream enters the mixing section at an angle of 2.5 
degrees. The splitter plate thickness has a base height of 0.5 mm at the trailing 
edge. Upstream of the straight sections for the two isolated flows shown in figure 
•i.l. contoured nozzle blocks were used to provide the supersonic flows with nearly 
uniform exit flow conditions. 

The mixing section height was 18 mm, and the overall length of the mixing section 
available for flowfield measurements was 500 mm. The width of the mixing section 
was 90 mm. and as a result, the mean flow development could be considered 
two-dimensional. I his was also verified in the experiment. The divergence angle of 
the lower and upper walls of the mixing section were adjusted in each experiment 
with two incoming supersonic flows, to account for boundary layer growth along 
these two surfaces and to effectively remove any streamwise pressure gradient. 

Single component LD\ measurements of turbulence intensities in the upstream 
flows taken 2 mm upstream of the splitter plate trailing edge indicated that 
incoming boundary layers were turbulent for all cases. These LDV measurements 
were also used to calculate the boundary layer, displacement, and momentum 
thicknesses of the two streams as they enter the mixing section. This makes the 
Goebel-Dutton experiments one of the more thoroughly documented benchmark 
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Figure 5.1: Schematic of Goebel-Dutton mixing layer experiment 


data sets available for compressible mixing layers. In the mixing region, a 
two-component LDV system was used to measure the axial and ti ansvei se 
velocities. In addition, a Schlieren system with a 20 ns pulse duration was used to 
obtain nearly instantaneous snapshots of the mixing layer. 

Goebel and Dutton examined seven cases using this experimental configuration. 
This work investigates their case 2 experiment. The operating conditions of the two 
streams in case 2 are provided in table 5.1. The two supersonic flows were matched 
in static pressure at the beginning of the mixing section (end of the splitter plate). 
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top flow bottom flow 


Mach Ye. 

1.91 

1.36 

1 '(»t/s) 

700 

399 

Tt(h') 

578 

295 

T(K) 

334 

215 


366 

293 

P(kPa) 

49 

49 

p[kgjm :i ) 

0.51 

0.79 

8 (mm ) 

2.9 

2.5 

8 m (mm ) 

0.90 

0.44 

0(mm) 

0.29 

0.21 


Table 5.1: Flow conditions for case 2 of the Coebel-Dutton experiments 


5.2 Two Dimensional Computational Modeling 

The development of the computational model began bv using the results of the 
Mach 1.91 and Mach 1.36 boundary' layer simulations discussed in chapter 4. 
Specifically, the wall-function solutions were examined, because the computational 
grids utilized in conjunction with the wall-function approach were constructed to 
enable a. continuous grid into the mixing region for use with the hybrid RANS-LES 
solver. The objective was to construct two HANS regions that would provide 
boundary layer quantities £, 8 ", and 0 that nearly matched those measured in the 
experiment, and shown in table 5.1. Because it would be virtually impossible to 
match all of three quantities exactly, the momentum thickness 0 was chosen as the 
key boundary layer parameter to match the computations with the experiment. The 
momentum thickness represents the mean momentum deficit entering the mixing 
section and is fundamental to the downstream mixing layer behavior. 
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Goebel-Dutton Ex periment Wall-Function Calculations 

^1 .91 (mm ) 

^i. 9 i {mm) 

Si.m{inw) 

S‘ l36 (mni) 

&\.m (mm) 

Table 5.2: Comparison of boundary layer quantities at splitter plate trailing edge 

While the comparison of calculated quantities to the Tucker theory in section 4.3 
was performed to verify the R.ANS method against well established correlations, the' 
objective here was to determine the axial length of plate needed for the Mach 1.91 
and Mach 1 .36 boundary layers to reach the same state as those' measured in the 
Goebel-Dutton experiment . Examining the Mach 1.91 boundary layer hist, table 
5.1 indicates that the momentum thickness for this stream was measured as 0.29 
rnw at the trailing edge of the splitter plate. For the wall function calculations 
discussed in section 4.3. the momentum thickness reached 0.29 mm at a Reynolds 
number of 3,540.000 (see figure 4.22), corresponding to an axial position of 198 mm 
from the leading edge of the plate. The other boundary layer quantities obtained 
from the calculations are also in close agreement with the Goebel-Dutton 
measurements, as shown in table 5.2. 

The next step was to construct a new computational grid. The axial domain was 
shortened to 198 mm while retaining the same number (141) of axial grid points. 
The grid stretching was modified to accommodate the shorter domain while 


2.9 

3.2 

0.90 

0.82 

0.29 

0.29 

2.5 

2.3 

0.44 

0.44 

0.21 

0.21 
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maintaining the initial and terminal grid s pacings. The vertical domain was reduced 
to 2T 70 mm to exactly match the height of the Mach 1.91 stream in the experiment. 
1 he vertical domain o( the original grid reached 22.75 mm at the 94th grid point in 
the vertical direction, so the first 94 points from the original grid were used in the 
modified grid. Calculations obtained with this 1 11 x 94 point. 198 mm by 22.75 mm 
grid provided boundary layer quantities identical to that of the original 141 x 141 
grid, further validating the grid independent characterist ics of the I?ANS method. 

A similar procedure was used for the Mach 1.26 boundary layer. Table 5.1 indicates 
that the momentum thickness for this stream was 0.21 mm at the splitter trailing 
edge. Examination of the wall-function solution obtained for the Mach 1.26 
boundary layer in section 4.2 revealed that the momentum thickness became 0.21 
mm at an axial position of approximately 120 mm, corresponding to a plate 
Reynolds number of 2.690,000. As was the case for the Mach 1.91 boundary layer, 
all three measures of t he boundary layer development were in close agreement with 
the experimental measurements for the Mach 1.26 case, as shown in table 5.2. 

A modified grid was constructed using 141 axial and 94 vertical points for this Mach 
1.26 case, corresponding to a physical domain of 120 mm by 22.75 mm. 

C alculations with this modified grid provided a solution identical to that obtained 
with the original 141 x 141 grid. 

The last step before constructing the entire RANS-LES computational grid was to 
use grid points from the modified wall-function grids just discussed, extending from 
81-141 in the axial domain and fix the solutions at the 81st axial station as the 
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inflow, in orclor to reduce the grid requirements of the final hybrid RANS-LRS giid. 
Boundary layer calculations subsequently obtained with these two shortened axial 
grids using 61 axial points arid 94 vertical points again returned the same boundary 
layer quantities and velocity profiles as the other solut ions. These were the RANS 
grids that were used to join with the LBS mixing region computational domain. 

As will be discussed in section 5.3, three different axial grid spacings were examined 
in the initial two dimensional hybrid calculations. In all cases, however, the vertical 
spacing from the two RANS regions was continued throughout the entire LES 
region. With the tightest vertical spacing of the wall-function boundary layer 
solutions set to 0.05 mm, 10 grid spacings were used vertically across the splitter 
base. As a result, all of the hybrid grids used 197 vertical points in the mixing 
region. 

For all of the hybrid calculations discussed in this chapter and in chapter 6. fixed 
inflow boundary conditions were used at the RANS inflows. The fixed inflow for the 
Mach 1.91 upper stream was placed at an axial position 67 mm upstream of the 
splitter plate trailing edge and the Mach 1.36 lower stream inflow was placed 42 mm 
upstream of splitter tip. Although information regarding the temperatuie and heat 
transfer characteristics of the splitter plate walls were not available from the 
experiment, these surfaces were set as adiabatic no-slip boundaries. At the outflow 
of the mixing section, corresponding to an axial distance of 300 mm from the 
trailing edge of the splitter plate, an extrapolation boundary condition was used, 
which is appropriate for the mixing supersonic flow exiting the axial domain. The 
top and bottom walls of the mixing section were approximated as slip walls, and no 
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attempts were made to simulate boundary layers developing on these two surface's. 
The Goebel-Dutton experimental configuration was specifically designed with 
adjustable divergence angles for these walls to account for boundary layer growth, 
and to thereby provide experimental data that may be directly compared to 
calculations which do not include the mixing section boundary layers. 

The combined RA.WS-LES calculations were obtained by marching in time with the 
Gottlieb- Turkel scheme for the entire computational domain. Although the R A X S 
regions did not change after reaching convergence, calculations were still performed 
in the RANS regions, in case any large perturbations from the LES region travelled 
upstream. Because of the zero-pressure gradient nature of these boundary layers, 
however, no large upstream fluct uations were noted in any of the hybrid 
calculations. No subgrid scale model was used in any of the two dimensional 
calculations discussed in this chapter. As a result, the use of a. turbulence model 
terminated at the end of the RANS regions. For the three dimensional calculations 
of the same configuration discussed in chapter 6. the Smagorinsky subgrid scale 
model was employed. 

For the two-dimensional cases discussed in this chapter and the three-dimensional 
cases discussed in chapter 6, a series of flowfield contours is shown in the mixing 
section to illustrate the mixing layer development and to help compare qualitative 
features of t he different modeling approaches. Instantaneous density and entropy 
contours are presented to show images of the mixing layer at a snapshot in time. 

For each of these quantities, an image is shown for the first one-third of the mixing 
duct to provide details of the initial mixing layer development. An image of the 
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entire computational domain is also shown, extending from the inflow of the two 
isolated HANS streams to the exit of the mixing section in the LES region. 

Density contours are useful for visualizing the flow characteristics in the mixing 
region. In addition, they provide a computational analogy to the Schlieren 
photographs that are used to illustrate mixing layer development in experiments. 
While the strongest gradients of the density will he observed to be in the developing 
shear layer, gradients are also observed in the regions above and below the mixing 
layer. These are the result of Mach waves generated by the unsteady mixing layer 
and their interaction with the two walls of the confined mixing section. Such waves 
were also evident in the Schlieren photographs of the Goebel-Dutton experiments. 

A Schlieren image taken of the first 250 rnw of the mixing section in the 
Goebel-Dutton case 2 experiments is shown in figure 5.2. 

In order to isolate the development of the mixing layer, instantaneous entropy 
contours are also used in the following discussions. Because the Mach waves 
occurring between the developing mixing layer and the two mixing section walls are 
of relatively weak strength, the entropy gradients in these regions are quite small. 

As a results, entropv contours enable emphasis to be placed on the mixing lavei 
development more clearly than is possible with the density contours. The entropy 
function of the PLOT3D program [109] is used here, which calculates the entropy in 
an ideal gas with constant specific heats as: 
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Y\ hero P r , j and y,< j are reference quantities used to nondimensionalize the entire 
flow domain. 

Time-averaged axial velocity and turbulent intensity profiles obtained from the 
calculations are compared with experimental measurements. The procedure used to 
obtain these time averaged quantities is discussed in appendix ( '. 

5.3 Axial Grid Density Effects 

Three computational grids are considered initially to examine the effects of axial 
grid resolution. Because the grid resolution in the vertical direction exceeded that of 
the axial spacing for even the finest grid examined here, no variation in the vertical 
spacing was considered and all of the grids had 197 vertical points in the mixing 
section. The three grids examined have 200, 400, and 800 axial points respectively 
and all are similar to the domain shown in figure 5.3, which corresponds to the 
coarsest grid used (200 axial points). The grid detail extending from the end of the 
RANS regions into the initial portion of the LES region is shown in figure 5.3(a). 
The entire computational domain is shown in figure 5.3(b). For clarity, only every 
third grid point is shown in both the axial and vertical directions in figure 5.3(b). 
Showing even' grid point for even the coarsest 200 point grid would obscure the 
depiction of the grid topology. 

Figure 5.4 shows the detail near the trailing edge of the splitter for all three grids. 

In every case, 10 equal grid spacings are used vertically across the 0.5 mm splitter 
plate trailing edge. This grid spacing (Ay = 0.05 mm) matches the initial vertical 
spacing of the wall boundary layers in the RANS regions. Axially, all three grids are 
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mixing region size 

200 axial pts. 

400 axial pts. 

800 axial | 

stretching factor 

1.0215 

1.0083 

1 .0029 

A.ri(/?m») 

.10 

.10 

.10 

A .rtast(i’i^i) 

6.52 

2.60 

1.06 


65.2 

26.0 

10.6 


Table 5.3: Comparison of axial grid spacitigs 

packed to the splitter trailing edge such that the initial A.r spacing is twice that of 
the finest A// spacing. The only difference between the three grids is the axial 
stretching factor, which is fastest for t he 200 point grid and slowest for the 800 
point grid. In all cases a geometric stretching factor was used. Table 5.3 provides a 
comparison of the axial spaeings for the three grids. 

The first, two dimensional simulation investigating grid density effects was for the 
coarsest grid using 200 axial points. Figures 5.5 and 5.6 provide instantaneous 
contours of the density and entropy for this case. In each of these contours, a vortex 
shedding pattern originates from the trailing edge of the splitter plate due to t he 
separation of the two flows leaving the wall bounded RAMS regions and entering the 
LES mixing section. The vortex shedding quickly dissipates and the flow appears to 
be laminar until an axial position corresponding to nearly one fourth of the overall 
duct length. At this position, an instability forms and the flow initiates transition to 
an unsteady turbulent pattern. The Reynolds number at the instability location is 
900,000 using an average value of the viscosity and density from the two streams, 
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tin' difference* in velocity of the two streams, and the axial position of approximately 
80mm. 


The t urbulent-like behavior quickly dissipates again, however, and the flow returns 
to a laminar state. For this case, the comparatively coarse grid enables the 
truncation error of the Gottlieb- Turkel scheme, shown previously in equation (3.24). 
to damp any oscillations without the use of any turbulence model or artificial 
dissipation in the mixing region. Because the turbulent behavior for this case was 
very limited, no turbulent averaging was done for this case. One final observation 
from this case was the Mach waves originating from the trailing edge of the splitter. 
These waves, in turn, reflect off of the mixing section walls and back onto the 
splitter. A qualitatively similar pattern was observed in the Goebel-Dutton 
Schlieren photographs as shown in figure 5.2. 

The second two dimensional simulation was for the computational grid using 400 
axial points. The characteristics of this flow were substantially different from those 
of the 200 grid point case. Instantaneous density and entropy contours are provided 
in figures 5.7 and 5.8. A stronger vortex shedding is evident for this case, although 
the vortex strength gradually dissipates back to a. laminar state approximately 40 
mm downstream of the splitter trailing edge. An instability again forms at an axial 
position of 80 mm which in turn slightly dissipates before resuming growth at a 
position 180 mm downstream of the splitter plate trailing edge. A turbulent pattern 
then grows from this location to the exit at x =300 mm. 
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The 400 point axial grid examined consisted of smaller axial grid steps than the 200 
axial point grid and had a significant ly reduced stretching factor relative to the 
coarse grid. As a result, the truncation error term in the Gottlieb- Turkel scheme 
that effectively smooths discontinuities is substantially reduced for t he 400 axial 
point grid, and the capability to resolve unsteady flow behavior is improved. 

The third computational grid investigated in this section was for the 800 axial point 
grid. A substantially different flow behavior is also observed for this case compared 
to the solutions obtained with 200 and 400 axial points. The density and entropy 
contours in figures 5.9 and 5.10 again indicate a vortex shedding pattern that 
originates from the trailing edge of the splitter plate, but unlike the other two cases, 
the solution does not return to a laminar state before transitioning over to a 
turbulent-like pattern. The very tight axial spacing for this case is sufficient to 
minimize the truncation error damping effects on the unsteady flow development. 
Interestingly, the transition from the organized vortex structure to a more random 
turbulent structure occurs at nearly the same location as the transition for the other 
two computational grids, although the flow behavior both upstream and 
downstream of this point is substantially different. 

A somewhat organized coherent structure may be observed from the contours of 
density and entropy for the 800 axial point case shown in figures 5.9 and 5.10. and 
near the end of the mixing section for the 400 axial point case shown in figures 5.7 
and 5.8. These structures are similar in form to those of the well known Brown and 
Roshko [16,86] mixing layer studies. The Brown and Roshko investigations of 
incompressible turbulent mixing layers indicated that the turbulent mixing layer 
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development was characterized by large coherent structures. As the Reynolds 
number in their experiments was increased, a fine scab' turbulence contained within 
the larger scale structures was evident . However, the mean flow properties were 
found to be the same regardless of the Reynolds number, indicating that the mean 
flow characteristics are dominated bv the large scale structure. 

The convective Mach number parameter was developed by Bogdanoff [12] and 
Papamoscshou and Roshko [76] for use in characterizing the compressibility 
characteristics of high speed turbulent mixing layers. For a planar mixing layer with 
equal specific heat ratios, the convective Mach number is defined: 
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For the Goebel-Dutton case 2 experiment investigated here, the convective Mach 
number was 0.455. Clemens and Mungal [21] conducted experiments of planar 
turbulent mixing layers with convective Mach number from 0.28 to 0.79. They 
found that as the convective Mach number is increased, the characteristics of the 
mixing layer changed from an organized two-dimensional Brown- Roshko structure 
to a three-dimensional structure with less evidence of large scale organization. 


Direct simulations of compressible shear layers conducted by Sandham and 
Reynolds [87] also indicted that three-dimensional instability modes become 
dominant at convective Mach numbers greater than 0.6. Their simulations also 
indicated reduced mixing layer growth rate with increasing convective Mach 
number. 
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The mean axial velocity profiles obtained from the two-dimensional solutions using 
400 and 800 axial grid points are compared to the measurements of Goebel and 
Dutton in figure 5.11. Because the calculations discussed in this chapter are only 
two-dimensional, and do not use a true LES approach in the mixing region, only 
qualitative comparisons with the data will be emphasized here. 4 he comparisons 
are made at four axial stations in the mixing layer, x = 50. 100, 150, and 200 mm. 
with x = 0 representing the beginning of the mixing section. As mentioned 
previously, x = 0 is also the axial position at which the switch from the BANS 
regions to the the LES region occurs. Comparisons of the two turbulence intensities 
u rms and v rma are shown in figures 5.12 and 5.13 respectively. 

The mean axial velocities are plotted versus vertical position normalized by the 
mixing section height, 48 mm. in order to give an indication of the shear layer 
spreading through the mixing section. In addition, the vertical positions were 
adjusted so that y/H = 0 represented the location where the local axial velocity was 
the mean of the two freestream velocities. The adjustment was also made in the 
experimental data as reported by Goebel and Dutton [36]. The turbulence 
intensities are plotted versus vertical position normalized by the local shear layer 
width, defined bv Goebel and Dutton as the distance between vertical positions 
where U = Ui - 0.1 AU and U = U 2 + 0.1 ML The two velocities f q and U 2 are the 
upper and lower stream velocities, respectively. 

The computed axial velocity contours indicate a much larger wake region than 
found in the experiments, and a greater discrepancy is observed for the 400 axial 
point case. At the x = 200 mm axial station, both solutions have returned to agree 
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moderately well with the data. The comparisons of st Teamwise and transverse 
turbulence intensities in figures 5.12 and 5.13 indicate substantial differences 
between the solutions and data. Both solutions, and especially the 400 axial point 
solution, show low levels of u rms and v rms at the .?■ = 50 nrm station. For the 400 
point case, this corresponds to the nearly laminar like state indicated by the 
contours shown in figures 5.7 and 5.8. Further downstream the intensities of t he two 
solutions grow substantially and by the x = 200 mm station, the transverse 
turbulence intensities from the solutions is substantially larger than those from the 
experiment . The profiles of ti rms obtained from the calculations demonstrates a 
double peak, particularly at ,v — 100 mm. This appears to be an indication that the 
velocity fluctuations are influenced more by the organized vortical structure early in 
the mixing sect ion than by the formation of turbulence. 

Despite the significant differences among the solutions obtained with the three 
computational grids discussed here, the vortex shedding from the trailing edge of 
the splitter and the appearance of waves originating from the splitter plate and 
reflecting off the mixing section walls were common to the three solutions. In the 
next section, the influence of the splitter base region and the mixing section walls 
are investigated. A final note from this section is that while grid refinement in 
HANS calculations is performed to achieve grid independence, and major changes in 
flow solutions when using two different grid resolutions usually indicates a modeling 
problem, this is not the case in LES computations. Grid refinement in LES enables 
more and more of the unsteady turbulent spectrum to be resolved, until the limit of 
DNS is reached. The purpose of a subgrid scale model in LES computations is to 
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Figure 5.2: Schlieren photograph of the Goebel-Dutton mixing layer (from Ref. :?1. 
used with permission) 


effectively represent the turbulent motion that is too small to be resolved by the 
computational grid. 
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(b) Entire computational domain (every third grid point in each direction shown) 


Figure 5.3: Computational grid for the 200 axial grid point case 
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(a) 200 axial points in mixing section 



(b) 400 axial points in mixing section 



(c) 800 axial points in mixing section 


Figure 5.4: Comparison of computational grids near the splitter plate trailing edg 
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(a) Beginning of mixing section 



(b) Entire mixing section 


Figure 5.5: Instantaneous density contours for the 200 axial grid point case 


NASA/TM — 2001-21081 1 


118 



(b) Entire mixing section 


Figure 5.6: Instantaneous entropy contours for the 200 axial grid point 
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case 


(a) Beginning of mixing section 



(b) Entire mixing section 


Figure 5.7: Instantaneous density contours for the 400 axial grid point case 
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(a) Beginning of mixing section 



Figure 5.8: Instantaneous entropy contours for the 400 axial grid point case 
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(b) Entire mixing section 


Figure 5.9: Instantaneous density contours for the baseline 800 axial grid point case 
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(a) Beginning of mixing section 



(b) Entire mixing section 


Figure 5.10: Instantaneous entropy contours for the baseline 800 axial grid point case 
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Figure 5.11: Time-averaged axial velocity profiles for 21) hybrid calculations investi- 
gating axial grid density effects 
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Figure 5.12: Profiles of n rm& for 2D hybrid calculations investigating axial grid density 
effects 
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Figure 5.12: Concluded. 
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Figure 5.13: Profiles of v rms for 2D hybrid calculations investigating axial grid density 
effects 
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Figure 5.13: Concluded. 
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5.4 Boundary Condition Effects 


In this sect ion. results obtained for the baseline case with 800 axial points are 
compared to two additional solutions obtained with 800 axial points. The first of 
these considers a modified splitter geometry in which t he splitter trailing edge is 
reduced to a sharp tip. and as a result, the flow separation and vortex shedding is 
removed. The second additional case is obtained for a modified mixing section in 
which the mixing section walls are moved very far vertically from t he mixing layer. 
This is done so that any waves originating from t he trailing edge of the splitter 
plate or the mixing layer do not have the opportunity to reflect off' the mixing 
section walls and back to the mixing layer. 

The overall grid structure for the case with a sharp trailing edge is identical to that 
of the baseline 800 axial grid point case, except for the treatment of the splitter 
plate trailing edge. A comparison of the grid detail around this region for the 
baseline geometry and the current case with a sharp trailing edge for the splitter is 
shown in figure 5.14. In the baseline geometry, 10 grid spacings are used in the axial 
direction at the base of the splitter, equally spaced at 0.05 mm. to resolve the flow 
region just downstream of the 0.5 mm thick splitter trailing edge. The modified 
geometry shown in figure 5.14(b), removes all but one grid spacing, such that the 
confined flows from the RANS regions will meet directly at the beginning of the 
LES region. The grid stretching in the vertical dimension was performed in the 
same manner as that used for the baseline 800 point grid, and reduced the vertical 
domain from 197 points to 188. The placement of the axial grid points was identical 
to that of the baseline grid. 
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In the Goebel-Dut ton experiment and in all of the computations discussed in the 
previous section, t lie height of the mixing section was 48 nun. I lie density contours 
shown in the previous section for the cases with 400 and 800 axial grid points, as 
well as Schlieren photographs taken in the experiment , revealed a series of Mach 
waves that reflected off of the top and bottom walls of the mixing sect ion and back 
onto the mixing layer. The strongest of these waves originated from the trailing 
edge of the splitter plate. 

To determine if these wave reflections influenced the instability formation in the 
shear layer where the simulations appeared to become turbulent, at j = 80 mm for 
both the 400 and 800 grid point cases, a modified grid was generated which moved 
each wall far from the mixing layer such the effective mixing section height became 
900 mm. This extreme spacing resulted in all waves generated from the splitter 
trailing edge or mixing layer to pass out of the outflow boundary at x — 300/??//? 
without the opportunity to reflect back onto the mixing layer. This computational 
grid, which had 800 axial points positioned the same as for the baseline 800 point 
case, is shown in figure 5.15. For clarity, only every fourth point is shown in each 
direction. In addition, the 197 points forming the vertical domain in the baseline 
grid were also used in this modified grid. The extra points needed to extend the 
vertical dimension to a total of 900 mm were added on to each side of the mixing 
section, such that the total number of points in the vertical dimension for this 
modified grid was increased to 347. 

Figures 5.16 and 5.17 provide instantaneous density and entropy contours for the 
case with a sharp trailing edge. As expected, the vortex shedding evident in the 
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baseline rase was removed in the current cast 1 wit h the sharp tip. The lack of a 
separation region results in an initially laminar mixing layer through the beginning 
of the mixing section. Interestingly, the laminar flow begins to transition to a more 
turbulent st ructure at nearly the same position observed for the baseline case, at 
approximately .r — 80 mm. The structures remain relatively small until 
approximately x = 150 mm where large scale t urbulence forms. These structures 
are more similar to the Brown-Roshko organized structures than were those of the 
baseline case with 800 axial points. Because the flows leaving the wall bounded 
regions experience a less rapid geometry change al the beginning of the mixing 
section for the sharp tip case, the strength of the initial waves off the splitter tip are 
reduced relative to the baseline case. 

For the second modified case with mixing section walls moved away from any region 
of influence on the mixing layer, the same set of instantaneous contours are provided 
in figures 5.18 and 5.19. Two thin lines are drawn on each of the contour plots to 
indicate where the mixing section walls were placed in the baseline case, however, 
the density contours for entire mixing section shown in figure 5.16(b) indicate that 
the flow domain extended beyond these thin lines. With the very large vertical 
domain indicated by the computational grid shown in figure 5.15, only a portion of 
the vertical domain is displayed in the contour plots. Several waves beyond those 
originating from the splitter tip are found to originate from the unsteady vortex 
shedding, but they do not reflect back to the mixing layer. 

As for the baseline 800 point case and the case with the sharp trailing edge, the 
transition location occurs at approximately x = 80 mm for this case, even without 


NASAyTM— 2001-210811 


132 



the influence of Mach wave reflection onto the shear layer. The downstream 
turbulent eddy structure is qualitatively similar to that for the baseline case, 
although the height of some of the structures is greater than t hat for the baseline 
case. This is likely the result of the lack of mixing section walls to confine the 
mixing laver for the modified case. One noticeable difference between the contours 
shown in figures 5.18 and 5.19 for the modified mixing section and those of the 
baseline case in figures 5.9 and 5.10 is that the mixing layer grows at a slightly 
upward angle without the presence of the mixing section walls. For the comparison 
of mean axial velocities and turbulence* intensities discussed next, an adjustment for 
the true mixing laver centerline will be employed, as was done in the previous 
section for the axial grid studies. 

In figure 5.20. mean axial velocity profiles obtained from the two-dimensional 
solutions investigating the different boundary treatments obtained with 800 axial 
grid points are compared to the Goebel-Dutton measurements. In addition, 
comparisons of the turbulence intensities u rms and v rms are made in figures 5.21 and 
5.22 respectively. As in the previous section, the primary objective here is to 
compare the different modeling approaches. Because these calculations were two 
dimensional only, and did not employ a subgrid scale model in the LES region, 
strict comparisons with the data are not emphasized. 

The mean axial velocity profiles at x = 50 mm indicate that the solution obtained 
with a sharp trailing edge had the smallest wake, which was expected because this 
case did not have the large separation and vortex shedding of the other two cases. 
Further down in the mixing section, all of the solutions are in reasonable agreement 
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with the data. The baseline case and t he solution obt ained with mixing section 
walls removed are in particularly close agreement with each other by the end of the 
mixing duct. 


All of the streamwise turbulence intensities are lower in peak magnitude than 
indicated by the Goebel-Dutton data at x = 50 nun. Further downstream in the 
duct, the magnitudes increase to be more in line with the data, and the two 
solutions obtained with the standard height for the splitter base exhibit a double 
peak. As discussed in the section 5.3, t his is believed to due to the influence of the 
upstream organized vortical behavior of the initial mixing layer. In figure 5.22. the 
solution obtained with a sharp trailing edge demonstrates significantly lower peak 
magnitudes in r rms at x = 50 mm and x = 100 nun. Further downstream in duct, 
all of the solutions indicate similar profiles of e rm5 , and all have significantly larger 
peaks than the experimental data. Lion et ah [59] and Inoue [45] also reported large 
over predict ions in the streamwise and transverse turbulence intensities for two 
dimensional mixing layer computations and attributed the discrepancies to the lack 
of a third computational direction. 

The results of this section indicate that the geometry of the splitter plate trailing 
edge has a significant effect on the initial shear layer formation. The unsteady 
vortex shedding is fundamental to the formulation of the hybrid RANS-LES method 
because although the mean flow properties of the incoming boundary layers are 
provided by the RANS regions, no turbulent oscillations are imposed that may 
initiate the instability of the mixing layer. Investigations of the splitter base region 
conducted by Clemens and Mungal [21] for an experimental configuration very 
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similar to that considered here also indicated an initial vorlex structure that 
transitions into turbulence. This transition, however, occurred significantly closer to 
the splitter plate trailing edge in the experiments than indicated by these two 
dimensional calculations. The three dimensional calculations discussed in the next 
chapter investigate this transition position furt her. 

Comparing the solution obtained with the standard placement, of the mixing section 
walls to the solution obtained with the walls moved vertically to prevent any Mach 
wave reflection back onto the shear layer, some minor differences were noted. 
However, there was little effect on the transition location from the organized vortex 
structure to turbulence. Because the Mach wave behavior observed in the baseline 
calculations was very similar to that indicated by the Schlieren photographs taken 
from the experiment, the standard placement of the mixing section walls will be 
investigated for the three dimensional calculations discussed in chapter 6. 
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(a) baseline grid 



(b) modified grid with sharp trailing edge for the splitter 


Figure 5.14: Comparison of computational grids using 800 axial points, near the 
split ter plate trailing edge 
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Figure 5.15: Computational grid for the 800 axial grid point case with the 
section walls removed (every fourth point shown in each direction) 
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(a) Beginning of mixing section 



(b) Entire mixing section 


Figure 5.16: Instantaneous density contours for the 800 axial grid point case with a 
sharp trailing edge for the splitter 
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(a) Beginning of* mixing sect ion 



(b) Entire mixing section 


Figure 5.17: Instantaneous ent ropy contours for the 800 axial grid point case with a 
sharp trailing edge for the splitter 
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(a) Beginning of mixing section 



(b) Entire mixing section length only, entire vertical domain is not shown 


Figure 5.19: Instantaneous entropy contours for the 800 axial grid point case with 
mixing section walls removed 
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(b) A — 100 mm 


Figure 5.20: Time-averaged axial velocity profiles for 2D hybrid calculations investi- 
gating boundary condition effects 
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Figure 5.21: Profiles of u rms for 2D hybrid calculations investigating boundary 
dition effects 
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Figure 5.21: Concluded. 
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Figure 5.22: 
dition effects 
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Figure 5.22: Concluded. 
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CHAPTER 6 


THREE DIMENSIONAL MIXING LAYER 
CALCULATIONS 


In the previous chapter, the two-dimensional calculations were 4 used to construct the 
initial computational model of the mixing layer and to examine preliminary effects 
of grid resolution and boundary condition treatment. To correctly investigate the 
capability of the hybrid method, however. LES calculations obtained in three spatial 
directions with the use of a subgrid scale model are required. I hese calculations are 
the focus of this chapter. The procedure used to extend the two-dimensional 
computational model described previously in section 5.2 to three dimensions is 
presented first in section 6.1. The results of three dimensional calculations obtained 
with the hybrid RANS-LES method are presented in section 6.2. 

6.1 Three Dimensional Computational Modeling 

The grid topologies and boundary conditions used for the three dimensional 
simulations were very similar to those used for the two-dimensional simulations 
discussed in chapter 5. The three computational grids with 200. 400, and 800 axial 
point grids described in chapter 5 were used to construct the three dimensional 
grids used here. To add the third computational direction in each case, the two 
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dimensional planar grid was copied to provide 11 points in the third (or z) direction. 
As a result, a side view of the three dimensional grids is represented by figure 5. A 
The grid spacing in the z direction was uniform and set equal to the axial spacing at 
the splitter trailing edge, or Ac = Ajq = 0.10 mm. Because of the very small 
number of grid points used in the c direction and the small physical space that is 
represented, only very small wave components in t his direction could be simulated, 
arid a periodic flow is assumed in this direction. 

The periodic boundary condition used in the c direction allows waves passing 
through one side of the ~ domain to enter the other side. Because t he 
Gottlieb- Turkel predictor corrector scheme uses a five point centered stencil in each 
direction, points along each boundary and one point interior to each boundary must 
be obtained when periodic boundary conditions are used. The solution vectors Q 
shown for the RANS and LES equation sets in equations 3.8 and 3.14 respectively 
are updated along the boundary corresponding to k — 1 in computational 
coordinates (Tj, k) as: 


Qay 1 — 3 

Qj,j,2 — fcmax— 2 

Similarly, along the other extreme boundary, corresponding to k 
solution vectors Q are updated using: 


( 6 . 1 ) 


kmax, the 


Qi, j./emax— 1 Q/J.3 

Q ijj.kmax — Qt,./\4 


( 6 . 2 ) 
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All of the other boundary conditions and the solution procedure are identical to 
that used for the two dimensional mixing layer calculations discussed in chapter 5. 
with the one exception being that the Smagorinsky subgrid scale model was used in 
these three dimensional simulations. As discussed in section ).2. no subgrid model 
was used in the two dimensional mixing layer calculations. 

The switch from the RANS regions to the LES region at the mixing plane 
(corresponding to a vertical plane drawn through the trailing edge of the splitter 
plate) is accomplished by changing the eddy viscosity used in the flow solver from 
the Cebeci-Smith turbulence model to the Smagorinsky subgrid scale model. As a 
result, the effect of the eddy viscosity changes from that of replacing all of the 
turbulent stresses in the RANS regions to that of only replacing the subgrid stresses 
in the LES regions. 

6.2 Three Dimensional Simulations 

The first, three dimensional simulation was obtained using the computational grid 
with 200 axial points. The Smagorinsky subgrid scale model used the coefficient 
C s = 0.1 and the standard expression for the subgrid length scale indicated 
previously in chapter 2 is repeated here: 


A = (AxAyAz)* (6.3) 

Figures 6.1 and 6.2 show instantaneous density and entropy contours at the middle 
plane in the ~ direction for these initial three dimensional simulations. With the 
very small domain in the 2 direction, the contours on all of the two-dimensional 


NASA/TM — 200 1-210811 


151 



x — if planes are very similar in appearance. The behavior of t he flow just 
downstream of the splitter trailing edge is very similar to that obtained for the two 
dimensional calculations shown in figures 5.5 and 5.6. The lack of adequate axial 
grid resolution results in a rapid dissipat ion of the initial vortex pattern, and is even 
more rapid in the three dimensional case due to the dissipative nature of the 
Smagorinskv subgrid scale model. Further out in the mixing section, the three 
dimensional calculations show no evidence of a secondary unsteadiness that was 
evident in the two-dimensional calculations, which is also due to the dissipation of 
the subgrid scale model. The shear layer appears to be effectively laminar 
downstream of the initial vortex region. Because no turbulent behavior was 
observed in the mixing section, no turbulent averaging was performed for this case. 

The three-dimensional grid using 400 axial points was utilized next, with 
substantially different flow development observed than that for the 200 axial point 
grid. As a result, three computations were performed in which the Smagorinskv 
constant and the subgrid model length scale were varied. The first two cases used 
the standard length scale expression given in equation (6.3) and investigated the 
two quoted extreme values for the Smagorinskv constant, C H ~ 0.10 and C s = 0.24. 
The third case set C s to 0.24 but used the modified length scale expression 
previously shown in chapter 2. and repeated here: 


A = 


(A.r) 2 + (A/yr -HAcf d 
3 


(6.4) 


For a computational grid with Ax — Ay — A z, equations (6.3) and (6.4) will return 
the same value for the Smagorinskv model length scale. The motivation for using 
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the modified expression was for the case of significant grid stretching, as occurs in 
both the axial and vertical directions away from the splitter plate trailing edge. 

Instantaneous contours of density and entropy at the middle plane in the ; direction 
are provided in figures 6.3 and 6.4 for the first case with the standard length scale 
expression and ( \ = 0.10. These contours indicate a fundamentally different flow 
structure than that observed for the two dimensional cases. In particular, the vortex 
shedding is observed to disintegrate into a random turbulent pattern much closer to 
the trailing edge. While the Sehlieren photographs taken in the Goebel-Dutton 
experiments did not probe the flow details near the splitter plate trailing edge, such 
details were examined by Clemens and Mungal [21] for a very similar experiment 
configuration and flow operating conditions. The Sehlieren image shown in figure 
6.5 indicates an initial flow structure similar to that observed in the calculation with 
an initial vortex shedding from the trailing edge of a splitter plate followed by a 
transition into turbulence. The turbulent structure in the experiment is also 
observed to be of primarily small scales, while the LES calculations, bv definition, 
only capture the large scale structures. While the length of the organized vortex 
structure in the calculations does not exactly match that of the Clemens-Mungal 
experiment, qualitatively the three dimensional calculations are in much closer 
agreement than the previous two dimensional results. In addition, the splitter plate 
thickness was 0.8 mm in the experiment of Clemens and Mungal, while that 
modeled in the calculations w^as 0.5 mm, which may be responsible for the precise 
differences between the experiment and calculations. 
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I he spreading of the turbulent shear layer is observed to be greater in this three 
dimensional case than in any of t he two dimensional cases. This may be observed in 
the instantaneous contours. A direct comparison of the density contours for the two 
and three dimensional cases obtained with 400 axial points is provided in figure (i.(i 
near the splitter tip. This comparison clearly shows the very different flow structure 
revealed by the three dimensional computations. Although the extent of the c 
domain is very small compared to the height of the mixing section, and cannot 
capture large structures in this direction, an unsteady mechanism in which 
disturbances may develop in all three directions and result in the rapid transition to 
turbulence is enabled by these three dimensional calculations. The two dimensional 
calculations, bv their very nature, do not allow for such three dimensional 
disturbances to develop. These results verify that LES calculations must be run in 
three dimensions in order to properly describe the initial turbulent flow structure. 

Time series snapshots of density contours and entropy contours immediately 
downstream of the splitter plate trailing edge are provided for this initial three 
dimensional case with 400 axial points in appendix D. The breakdown of the 
organized vortical structure originating from the splitter plate wake into a turbulent 
structure is demonstrated in these time series. A mentioned previously, the 
transitional behavior near the trailing edge of the splitter was not evident with the 
two dimensional calculations. The frequency of the vortex shedding was the same 
for the two and three dimensional cases. An analysis of this shedding frequency, also 
provided in appendix D. shows that a Strouhal number calculated using the splitter 
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base height and difference in velocities of the two streams is 0.225. which is very 
close to t he Strouhal number range characteristic of flow past cylinders, 0.20 - 0.21. 

Examining the other two solutions for the three dimensional approach with 400 
axial gr id points next, the instantaneous density and entropy contours for the case 
using the standard length scale and C s = 0.24 are provided in figures 6.7 and 6.8. 
The same contours for the case using the modified length scale expression and 
C\ = 0.24 are shown in figures 6.9 and 6.10. Although the instantaneous contours 
for these two cases and the initial 400 point case differ in exact structure in the 
snapshots shown, overall the turbulent structures are qualitatively the same for the 
three cases. In addition, the number of organized vortices before breakdown to 
turbulence may be observed to differ for these three cases at the particular instants 
in time shown, but in examining each of the solutions over development, in time, 
each of the solutions oscillated in having between 5 and 10 organized vortices. 

The mean axial velocity profiles obtained from the three dimensional hybrid 
RANS-LES calculations using 400 axial grid points in which the subgrid model 
parameters were varied are compared to the Goebel- Dutton data, in figure 6.11. A 
comparison of the two turbulence intensities is made in figures 6.12 and 6.13 
respectively. Examining the mean axial velocities first, all of the three dimensional 
solutions exhibit a larger wake at x = 50 mm than that of the experiment, but 
compared to the two dimensional solutions, the three dimensional wakes are 
significantly smaller. This improvement is the result of the more turbulent behavior 
for the three dimensional cases in the beginning of the mixing section, further 
downstream, the three dimensional solutions indicate reasonable agreement with the 
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data, although all of the solutions appear to mix more rapidly than indicated by the 
experimental data. 

The profiles of u rms in figure 6.12 generally indicate overpredictions from the 
calculations, which corresponds to the wider axial velocity profiles in figure 6.1 1. At 
x = 50 77 ?;??. a double peak in the calculated intensities is somewhat evident, 
although the effect is much less pronounced for these three dimensional calculations 
than for the two dimensional calculations in chapter 5. Further downstream, the 
three dimensional solution obtained with the modified length scale demonstrates 
generally lower levels of u rms than the other solutions. The modified length scale 
expression results in larger values of the subgrid model eddy viscosity, which in turn 
damps more of the small scale unsteadiness. Although increasing the Smagorinskv 
constant C s also tends to result in more damping, the effect of changing C s from 
0.10 to 0.24 does not seem to have as large an effect as the subgrid scale length 
expression. 

The computed profiles of r rma , shown in figure 6.13 are also overpredicted relative to 
the data, with the lowest levels of v rms predicted with the modified length scale and 
C 8 = 0.24. As mentioned previously. Liou et al. [59] and Inoue [45] also reported 
large overpredictions in the turbulence intensities for planar mixing layers. With the 
confinement of the current three dimensional calculations to a very small domain in 
the c direction, the inability to calculate large scale fluctuations in this direction 
may be responsible for the overpredictions of u rrns and v rms . In addition, the 
Schlieren photographs taken in the Goebel- Dutton experiment indicated a very fine 
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turbulent structure contained within the larger scale development , while the EES 
calculations inherently can resolve only the larger turbulent scales. 

A final three dimensional case was run using 800 axial points in the mixing section. 
The modified length scale expression and ( = 0.24 was used for the subgrid model. 
The finer grid used here reduced the permissible time step by nearly a factor of two 
relative to the 400 axial point cases. This reduced time step and doubling the 
number of grid point s in the mixing section would require a factor of four increase in 
the computer CPU time requirements to run this case to completion, relative to the 
cases with 400 axial points. Considering that each of the 400 axial point cases 
required 500 CPI hours oil a ( ray C90 computer, 2000 ( ray ( 00 hours would be 
required to complete the case with 800 axial grid points. As a result, this last three 
dimensional case was run long enough to allow the flow to fully develop in the 
mixing section, but not long enough to enable time averaging of the turbulent 
statistics. 

In contrast to the grid refinement studies performed for the two dimensional cases, 
the large scale turbulent development did not change significantly when increasing 
the number of axial points in the three dimensional computations from 400 to 800. 
The instantaneous density and entropy contours in figures 6.14 and 6.15 indicate a 
large scale turbulent structure which closely resembles those of the 400 axial grid 
point case. In particular, the breakdown of the organized vortex structure to 
turbulence is very similar to those previously shown for the 400 axial grid point 
case. Within the large scale structures, more fine scale turbulence is evident in the 
800 axial point case. This behavior corresponds directly with the philosophy of LES 
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in that as the computational grid is refined, smaller structures are able to be 
resolved and the role of the subgrid scale model is reduced. In the idealized limit of 
an infinitely fine grid, a direct numerical simulation is obtained. 
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(a) Beginning of mixing section 
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(b) Entire mixing section 


Figure 6.1: Instantaneous density contours for the 200 axial grid point case (3D) 
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(b) Entire mixing section 


Figure C.2: Instantaneous entropy contours for the 200 axial grid point case (3D) 
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(b) Entire mixing section 


Figure 6.3: Instantaneous density contours for the 400 axial point case using the 
standard turbulent length scale and C s = 0.10 
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(a) Beginning of mixing section 



(b) Entire mixing section 


Figure 6.1: Instantaneous entropy contours for the 400 axial point case using the 
standard t urbulent length scale and C s = 0.10 
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(a) Beginning of mixing section 



(b) Entire mixing section 


Figure 6.7: Instantaneous density contours for the 400 axial point case using the 
standard turbulent length scale and C s = 0.24 
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(b) Entire mixing section 


Figure 6.8: Instantaneous entropy contours for the 400 axial point case using the 
standard turbulent length scale and C a — 0.24 
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(b) Ent ire mixing section 


Figure 6.9: Instantaneous density contours for the 400 axial point case using the 
modified turbulent length scale and C s = 0.24 


N AS A/TM— 200 1-210811 


167 


(a) Beginning of mixing section 



(b) Entire mixing section 


Figure 6.10: Instantaneous entropy contours for the 400 axial point case using the 
modified turbulent length scale and C s = 0.24 
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Figure 6.11: Time-averaged axial velocity profiles for 3D hybrid calculations 
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(b) A" = 100 mm 


Figure 6.12: Profiles of u rrns for 3D hybrid calculations 
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Figure 6.12: Concluded. 
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Figure 6.13: Profiles of v rms for 3D hybrid calculations 
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Figure 6.13: Concluded. 
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(b) Entire mixing section 


Figure 6.14: Instantaneous density contours for the 800 axial point case using the 
modified turbulent length scale and C a = 0.24 
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(b) Entire mixing section 


Figure 6.15: Instantaneous entropy contours for the 800 axial point case using the 
modified turbulent length scale and C s = 0.24 
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CHAPTER 7 


CONCLUSIONS AND RECOMMENDATIONS 

The work described in this dissertation represents the initial efforts to develop and 
evaluate a hybrid RANS-LES method for compressible mixing layer simulations. 
Such mixing layers dominate the flows in exhaust, systems of commercial and 
military aircraft in current use and also those of hypersonic vehicles under 
development for future space transportation use. The hybrid method uses a RANS 
approach to provide the mean flow characteristics of the wall boundary layers 
entering the mixing layer and an LES approach for the mixing region. Although the 
RANS approach does not provide any unsteady turbulent information to the LES 
region, the mean flow boundary layer characteristics are provided. The hybrid 
method was developed for the analysis of nozzle and mixing layer configurations in 
which a dominant structural feature, such as the base region of a nozzle or splitter 
plate separating the upstream flows, will provide the dominant unsteady mechanism 
to drive the development of turbulence in the mixing layer. 

The hybrid method development was initiated by deriving a set of 
Reynolds-averaged Navier-Stokes (RANS) equations using density weighting in the 
averaging process, and a set of spatially-filtered large eddy simulation (LES) 
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equations which also used density ( Favre) weighting. The resulting similar form of 
the RANS and LES equation sets enabled both to be solved with a single solution 
scheme. In this dissertation, the Got tlieb-Turkel predictor-corrector scheme was 
employed. A transformation of the equations to generalized coordinates enabled 
flow calculations on stretched, non-Cartesian grids. The RANS equations were 
closed using the Cebeci-Smith algebraic turbulence model, with the option to 
employ the wall-function technique of Ota and Goldberg. The LES equations were 
closed using the Smagorinsky subgrid scale model. 

The Cebeci-Smith turbulence model, despite its relatively simple form, was 
demonstrated to provide accurate calculations of boundary layer flows that are free 
of adverse pressure gradients or separation regions. Further, the use of the 
Cebeci-Smith model in conjunction with the Ota-Goldberg wall function enabled 
calculations of supersonic wall boundary layers to nearly the same accuracy as that 
of the standard approach of integrating the Cebeci-Smith model through the viscous 
sublayer, while enabling a significantly larger vertical grid spacing near the wall. As 
a result, the wall function approach enabled a continuous computational grid to be 
used from the RANS to the LES regions, and the method thereby avoided the use of 
discontinuous grid zones that would have otherwise required an interpolation 
scheme between the two regions. In addition, the origins of the Cebeci-Smith RANS 
turbulence model and the Smagorinsky LES subgrid scale model are both in mixing 
length theory, and this similar form of the two models assisted in code 
implementation. As a result, the use of a more sophisticated turbulence model to 
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close the RANS equations was found to be unnecessary, provided the HANS regions 
are restricted to attached, zero pressure gradient wall boundary layer regions. 

While true LES calculations require computations in three spatial directions, two 
dimensional simulations of a benchmark mixing layer experiment were considered 
first to address the effects of axial grid resolution and boundary conditions. The 
parametric study of axial grid resolution indicated more realistic turbulent 
development with increasing axial grid density. For the coarsest grid examined, 
there was almost no evidence of turbulent flow development. For all of the case's 
examined, a vortex shedding was found to originate from the base region of a 
splitter plate separating the upstream wall bounded regions. For the finest two 
dimensional grid examined, the unsteady vortex pattern eventually transitioned to a 
turbulent structure. The location of this transition, however, was much further 
downstream than observed in the experiments. 

Additional two dimensional calculations were obtained to investigate the boundary 
treatment of the splitter plate trailing edge and of the mixing section walls. 
Calculations obtained for a case in which the finite thickness splitter base was 
changed to a sharp tip indicated that the vortex shedding was removed, but the 
development of turbulence downstream occurred at nearly the same position as for 
the case with vortex shedding induced by the splitter base geometry. Computations 
were also obtained for a modified geometry in which the mixing section walls were 
effectively removed, to determine if Mach waves reflecting off these walls in the 
baseline calculations affected the turbulent mixing layer development. These 
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calculations indicated that removing the* mixing section walls did not change the 
fundamental structure of the flow relative to the baseline case. 

Three dimensional calculations were obtained next for grids constructed by copying 
the two dimensional planar grids to locations in the third computational direction, 
normal to both the streamwise and transverse directions. Only a small domain was 
modeled in this third direction, and periodic boundary conditions were employed 
along the extreme boundaries. For the coarse three dimensional grid, again no 
turbulent flow development was observed. For the intermediate grid, the vortex 
shedding found previously in the two dimensional simulations was also observed in 
the three dimensional calculations. However, the organized vortical structure 
rapidly disintegrated into a significantly more realistic turbulent flow structure. 

This rapid transition to turbulent flow was nearly identical to that found in 
experimental investigations of a similar mixing layer configuration. Although the 
extent, of the third dimension in these calculations was very small, an unsteady 
mechanism by which disturbances could develop in all three directions and result in 
a rapid transition to turbulent flow was enabled by the three dimensional 
calculations. In contrast, a two dimensional approach, bv definition, does not. allow 
for such three dimensional disturbances to develop. The results of these calculations 
verified that LES simulations must be performed in three dimensions. 

Parametric studies of the subgrid model length scale and the Smagorinsky model 
coefficient were examined with this intermediate grid, but no significant differences 
were noted. Comparisons of time-averaged axial velocities and turbulence intensities 
from the calculations to experimental data indicated reasonable agreement, with the 
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solutions indicating somewhat higher levels of turbulent mixing. A ma jor source of 
discrepancy between the calculations and experiment is believed to be the lack of 
adequate grid resolution to resolve the small turbulent scales contained within the 
larger turbulent structures. Another source of discrepancy was the very small 
domain used in the third computational direction. 

Despite these limitations, the three dimensional calculations demonstrated the 
success of the hybrid method to capture the dominant characteristics of the mixing 
layer, and in particular, the rapid transition of the organized vortex structure to a 
turbulent mixing layer structure. It is expected that improvements in the fidelity of 
the solution scheme, and more importantly, improvements in computing power, will 
enable better predictions of the turbulent statistics, as will be discussed briefly. 

A final three dimensional calculation was investigated using a computational grid 
constructed from the most densely packed two dimensional grid. Because a 
prohibitively long run time would be required to complete this solution for turbulent 
statistics purposes, the calculations were run only long enough to allow the flow to 
fully develop in the mixing section. The large scale turbulent structures evident for 
this case were very similar to those for the intermediate three dimensional cases. 
More resolution of the finer turbulent scales contained with the larger structures 
was observed for the fine grid case, which is in line with the philosophy of LES to 
resolve finer scales as the grid density is increased. 

The effects of improved subgrid scale models on the quality of LES simulations is an 
issue of considerable debate. Research into advanced subgrid scale models has 
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yielded improvements in the accuracy of LES simulations obtained for low Reynolds 
number wall boundary layer flows. However, authors such as Spalart [96] and 
Fureby and Crinstein [29] offer the opinion that subgrid scale refinement will offer 
only small improvements for large eddy simulations of flows away from boundaries 
and without chemical reactions. Both aut hors further suggest that it may even be 
feasible to perform an LES simulation without an explicit subgrid scale model, 
provided the numerical scheme is sufficient to prevent unresolved wavenumbers from 
contaminating the solution and that the simulation resolves turbulent scales in the 
inertial subrange. 

Improvements to the numerical method may enable more accurate resolution of the 
smaller turbulent scales, which appeared to be the greatest limitation of the current 
method in simulations of the compressible mixing layer experiment. The class of 
numerical methods known as compact schemes, such as those presented by Lele [55]. 
have been shown to provide higher spatial accuracy and improved capability to 
resolve higher wave numbers for a given grid size than is possible using 
MacCormack-tvpe predictor-corrector schemes such as the Gottlieb- Turkel method. 

The continuing advances in computing speed and computer memory will also enable 
calculations of higher fidelity, arid eventually to the limit of direct numerical 
simulations (DNS) where all turbulent scales of importance are resolved. Estimates 
of Spalart [96]. however, suggest that improvements of several orders of magnitude 
in computer speed and memory will be required to perforin full LES or DNS 
calculations of realistic engineering configurations. His estimates project full LES 
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simulations requiring 1 0 1 1 grid points to be possible' in t Ire year 2045. and DNS 
simulations requiring 10 1 ' 1 grid points to be possible in the year 2080. 

In the interim. RANS and hybrid RANS-LES methods will each be the most 
appropriate computational tool for certain classes of turbulent flows. For wall 
bounded turbulent flows without massive separation regions. RANS methods will 
likely be the most appropriate choice for some time. Hybrid methods, such as that- 
developed in this work, will very possibly become the appropriate tool for flows with 
significant mixing regions or large scale separation zones. As a result, research to 
further develop both RANS and hybrid RANS-LES methods will be important to 
improve the capability to simulate complex turbulent flows. 
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APPENDIX A 


ACCURACY ANALYSIS OF THE GOTTLIEB-TURKEL 

SCHEME 


In this appendix, the temporal and spatial accuracy of the Gotti ieb-durke] scheme 
are investigated. 

The Gottlieb- Turkel scheme is illustrated using a one dimensional problem that is a 
model for the Navier-Stokes equations written in vector form. This model problem 
is given as: 


Oq 

dt 


df 

d.l 


(A.l) 


The predictor step is: 


( li = </« 
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The corrector step is: 
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The time step of the Gottlieb-Turkel predictor-corrector scheme, A. is related to the 


grid spacing, A,r, and the propagation speed, .4. through the CFL 
(( 'ourant-Friedrichs-Levvy ) number: 


At = CFL^f (A. 1) 

A 

Nelson [71] investigated the spatial and temporal accuracy of several 
predictor-corrector schemes, including the Gottlieb-Turkel scheme, investigated here 
using a linearization of the model problem shown in equation A.l. This procedure 
will also be used here. The linearization of equation (A.l) begins by assuming a 
constant propagation speed. A: 


f = Aq (A. 5) 

The equation for the predictor step (A. 2) can be substituted into the equation for 
the corrector (A. 3). which results in 
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The terms in this equation having a superscript (*). which represent the fluxes / 
after the predictor update, can be rewritten using /* = Aq*. Equation (A. 6) then 
becomes 
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The q* terms are replaced using the predictor expression (A. 2). 1 his enables the 
complete predictor-corrector update to be expressed as function of the original 
values of q and /: 


n +1 ii 

Qi = 


TTT — (~~.fi + 8 /,+ 1 — Ji+ 2 ) ~ 777 “T ( ~ ( h ~ ^ ( H - 1 + </:->) + 
12 A.r 1 2 A.r 


~ (x”” ) [*t — f fi C 8/<+ 1 — /i+'i) — S( T_/ ] + 8./ ( ./ i+l ) + ( ~Ji-2 + S / , _ 1 ./(')] 


Using / = .U/ and combining terms results in: 
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(A. 9) 


Next, the q* term is moved to the left side of the equation, and the q and / terms 
are expressed in terms of Taylor series expansions about the points q n and t /,-, 
respectively: 
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(A. 10) 


1 2 V A.r / W dx' 2 dx 4 

Reorganizing equation (A. 10) and dividing all of the terms by At results in: 


dq Of At d 2 q AC d 3 q Ax 4 cT f AAtAx 2 d 4 f AAt d 2 / 

at + “ 2 at 2 6 Ctt 3 + 30 fUr 5 IS <9.r 4 2 ch 2 

(A. 11) 


Finally, the first and fifth terms on the right side of equation ( A . 1 1 ) are removed by 
taking the second derivative of / = Aq. and then recognizing that 
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(A. 12) 


A i() 2 q _ A At d 2 f 
~2~df 2 ” 2 PC 

The resulting expression is: 


dq df A l 2 d 3 q A At Ax 2 d A f Ax A if f 

c)t rh () dP 18 dx 4 30 dr 5 1 ' 

The left side of equation (A. 13) is the original model problem given in equation 
(AT). As a result, the right side of equation (A. 13) is the truncation error resulting 
from the discretization of the model equation using the Gottlieb- Turkel scheme. 

The first term on the right side of equation (A. 13) indicates that temporal accuracy 
of the scheme is second order. The second truncation error term involves both At 
and Ax 2 , and we can use relation between the time step, grid size and CFL number 
shown in equation (A. 4) to rewrite this equation (A. 13) as: 


dq df _ At 2 (Pq Ax 3 iff Ax A (P£ 

dt + dx 6 dP ^ 18 dx 4 ^ 30 Ox* 


(A. 14) 


The Gottlieb-Turkel explicit scheme is only stable for CFL values smaller than one 
For CFL values on the order of one, the spatial accuracy of the scheme is third 
order. In practice, the maximum CFL number is usually set to a value of 0.5 or 
less. For stretched grids, the limiting time step is proportional to the smallest grid 
spacing arid then the effective local CFL number is much smaller in regions of the 
computational domain where the grid spacing is larger. This can be observed by 
considering equation (A. 4) for the case of variable grid spacing. Ax, but a constant 
time step At. As a result, the truncation error term in equation (A. 14) becomes 
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insignificant away from the region of tightest grid spacing. For such regions, the 
next truncation error term indicates fourth order spatial accuracy. In conclusion, 
the Gottlieb-Turkel scheme is strictly second order accurate in time and third order 
accurate in space, but in the case of highly stretched grids, the spatial accuracy is 
effectively fourth order. 
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APPENDIX B 


TRANSFORMATION METRICS 


The metric terms that are required to transform the RANS arid LES equations from 
physical space to computational space are derived in this appendix. I he procedure 
used in this work is the same as that presented in reference [42]. 

The equations derived in chapter 2 are transformed from physical space {x.y,z) to 
computational space (£. ?/.(,') using the relations: 

£ = £U’,y-~) 

1/ = ?/(.r,.i/,r) (B.l) 

C = £(•<••!/•-) 

The chain rule of partial differentiation allows the cartesian derivatives to be 
expressed as: 
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In equation (B.2), the terms s \ . ?/.. and C. are the metrics of the 

transformation. To compute tliese metric terms, the first ste|> is to compute the 
derivatives .r,.-, ,r c . y y, r y < . c £ , ~, r and z c . The stepsizes of these derivatives are 
equal in computational space and can be obtained using finite difference 
expressions. To be consistent with the doth lieb-Turkel scheme, which effectively has 
fourth order spatial accuracy and uses a five point stencil, these spatial derivatives 
are also calculated with a five point, fourth order accurate finite difference method. 
Derivatives in x: 
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Derivatives in z: 
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r< =_ ~ 12 

For grid points that are located at a boundary, a one sided fourth order differenced 
expression is used, and at grid points one point off of a boundary, a skewed fouith 
order difference expression is used. Both expressions are obt ained in the same 
manner as the central differenced expressions in equations (B.3 - B.5) through use of 
the appropriate Taylor s series expansions. As an example, the ■>' c term at a 
boundary corresponding to (2 = 1) i s: 

— 25j’j 4~ 4(\r ;+ ] — Mxi'j +2 T 1 ~ J-Q-h ( B 6 ) 

Similarly, the x ( term at (j = 2) is 


— 'ixj-i — lO.r, + 18.t_/ + 2 — 6x^+3 + 't'i +3 


(BA 


Once all of these spatial derivatives are obtained, the Jacobian of the 
transformation. .7, is calculated: 


x t — x v {Vi-'c + x c (y«~» y < ?•*$) 
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Using this expression for the Jacobian, the metric terms are: 
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APPENDIX C 


TURBULENT TIME AVERAGING PROCEDURE 


The procedure used to obtain the time-averaged axial velocities and the two 
turbulent, intensities for the LES regions is detailed in this appendix. 

As each of the two dimensional calculations discussed in chapter 5 and the three 
dimensional calculations discussed in chapter 6 progressed in time, quantities were 
accumulated from the instantaneous velocity fields at each time step. These 
quantities were then used to calculate the mean axial velocity u and the two 
turbulence intensities u rms and v rms . The mean velocity field is calculated using: 


1 Jo 


(c.i; 


While the two turbulence intensities u rms and v rrns are obtained from: 


and 



(0.2) 



(0.3) 


NASA/TM — 200 1-210811 


195 



The two quantities t/' and v f are not known until 77 is determined, which in turn 
requires completion of the time averaging period T. As a result, an alternative to 
equations ((’.2) and ((M) is used to find the turbulence intensities -u rmB and r rms 
Using it = 77 + u f and r = r + rh equations (C.2) and (CM) can be rewritten as: 





For all of the two and three dimensional mixing layer simulations, a constant time 
step was used so that equations ((hi), ( ( 1 .4 ) . and ((M) can be obtained through a 
straightforward summation procedure. This procedure consists of storing 
summations of «. i\ t/ 2 , and v 2 at each grid point from each time step. Equations 
(CM), (CM), and (C.5) then can be replaced with equations (C.6), (C\8), and (C.9). 

1 A ' 

« = T7 5Z U (C.6) 

^ n = 1 

where V is the total number of time steps, corresponding to the total time interval 
7\ which are related through the time step size At: 


T = A At 
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In all of the mixing layer simulations, the initial iterations used to start the flow 
development were run using a CFL number of 0.3 wit h the Gott lieb-Turkel scheme 
As the unsteady flow developed, however, the minimum time step fluctuated in the 
flow, so that the actual time step given by equation 3.31 also fluctuated, even with 
the use of a constant CFL number. The smallest actual time step observed during 
these iterations needed to allow the flow to fully develop was monitored, and then 
this fixed time step was imposed for the iterations in which the turbulent statistics 
were accumulated. 


Each of the two-dimensional calculations were run for approximately three average 
flow'-through periods, once the init ial flowfield in the entire domain was established. 
This number of flow-through periods wras sufficient to allow the mean velocity and 
turbulence statistics to reach converged levels. Due to the larger turbulent 
unsteadiness that was found in the three dimensional calculations, these calculations 
were run for approximately four average flow-through periods, once the flowfield in 
the entire computational domain was fully developed. 
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APPENDIX D 


TIME SERIES OF INITIAL MIXING LAYER 
DEVELOPMENT 


In this appendix, a time series of density and entropy contours are provided to 
illustrate the initial mixing layer development for a three dimensional hybrid 
method calculation of the Goebel-Dutton experiment. The computation discussed 
here used the standard length scale expression and C = 0.1 with the Smagorinskv 
subgrid scale model. As discussed in chapter 5, the density contours provide a close 
analogy to the Schlieren photographic technique used in experiments to illustrate 
mixing layer development. In addition, entropy contours enable emphasis to be 
placed on the shear layer alone, and do not show the Mach wave contours between 
the shear layer and the mixing section walls that are evident in the density contours. 
These entropy contours are generated using the entropy function of the PLOT3D 
program. In reference [109] the entropy function is defined as: 

S-S„, =CJn{-^) +Cpln (^) (D.l) 

Figures D.l and D.2 show sixteen snapshots of the density and entropy contours, 
respectively. Each successive snapshot represents a march forward in time of 
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7.r 1 () —<> .seconds, which corresponds to 250 time steps obtained with the 
(iottlieb-Turkol predictor-corrector scheme. This time interval was chosen such that 
in every other snapshot, a new vortex is shed from the trailing edge of the splitter. 
The initial organized vortex pattern, which very similar to the Karman vortex street 
characteristic of the separated flow past a cylinder, extends approximately the 
length of six vortices before the mixing layer disintegrates into a more turbulent 
pattern. 

With the shedding of new vortices occurring approximately every 500 time steps 
with the numerical scheme, or 7.4.rl0“ 6 seconds, a Strouhal number may be 
calculated for further analogy with the separated flow past a cylinder, where in this 
case of the separated flow past the splitter plate: 


St = 


JH 

AU 


(D.2) 


Using the splitter plate base height of 0.5 mm and the difference in velocities of the 
two incoming streams equal to 300 mj . s, the calculated Strouhal number is 0.225, 
which is very close to the Strouhal numbers found for flow past cylinders [54,74] of 
0.20 - 0 . 21 . 
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Figure D.2: Concluded 
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three dimensions. Comparisons of time-averaged axial velocities and turbulence intensities indicated reasonable agreement with experimental data. 
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